Finite-element solution to multidimensional multisource electromagnetic problems in the frequency domain using non-conforming meshes

Finite-element solution to multidimensional multisource electromagnetic problems in the frequency... Summary We propose an approach to solving multisource induction logging problems in multidimensional media. According to the type of induction logging tools, the measurements are performed in the frequency range of 10 kHz to 14 MHz, transmitter–receiver offsets vary in the range of 0.5–8 m or more, and the trajectory length is up to 1 km. For calculating the total field, the primary-secondary field approach is used. The secondary field is calculated with the use of the finite-element method (FEM), irregular non-conforming meshes with local refinements and a direct solver. The approach to constructing basis functions with the continuous tangential components (from Hcurl(Ω)) on the non-conforming meshes from the standard shape vector functions is developed. On the basis of this method, the algorithm of generating global matrices and a vector of the finite-element equation system is proposed. We also propose the method of grouping the logging tool positions, which makes it possible to significantly increase the computational effectiveness. This is achieved due to the compromise between the possibility of using the 1-D background medium, which is very similar to the investigated multidimensional medium for a small group, and the decrease in the number of the finite-element matrix factorizations with the increasing number of tool positions in one group. For calculating the primary field, we propose the method based on the use of FEM. This method is highly effective when the 1-D field is required to be calculated at a great number of points. The use of this method significantly increases the effectiveness of the primary-secondary field approach. The proposed approach makes it possible to perform modelling both in the 2.5-D case (i.e. without taking into account a borehole and/or invasion zone effect) and the 3-D case (i.e. for models with a borehole and invasion zone). The accuracy of numerical results obtained with the use of the proposed approach is compared with the one obtained by other codes for 1-D and 3-D anisotropic models. The results of this comparison lend support to the validity of our code. We also present the numerical results proving greater effectiveness of the finite-element approach proposed for calculating the 1-D field in comparison with the known codes implementing the semi-analytical methods for the case in which the field is calculated at a large number of points. Additionally, we present the numerical results which confirm the accuracy advantages of the automatic choice of a background medium for calculating the 1-D field as well as the results of 2.5-D modelling for a geoelectrical model with anisotropic layers, a fault and long tool-movement trajectory with the varying dip angle. Electrical properties, Downhole methods, Numerical approximations and analysis, Numerical modelling 1 INTRODUCTION For processing the electromagnetic soundings data, the multidimensional modelling is increasingly used (e.g. Oldenburg et al.2013; Kamm & Pedersen 2014; Kordy et al.2016b). The inversions based on the 2-D or 3-D modelling make it possible to recover the structure of the investigated medium with greater adequacy (e.g. Persova et al.2013; Schwarzbach & Haber 2013; Haber & Schwarzbach 2014; Usui 2015; Persova et al.2016, etc.). However, these inversions require high computational costs. Therefore in practice, until recently, quite simple and low-cost approaches, though not always adequate, have often been used. They are based on calculating the electromagnetic fields in the 1-D medium or other approaches, which make use of simplified mathematical models of the electromagnetic field (among others Macnae 2015; Hunkeler et al.2016). In some cases, these approaches allow us to obtain acceptable results (e.g. Hennessy & Macnae 2010; Hunkeler et al.2016). Nevertheless, the researchers often note that the simplified and 1-D approaches not rarely give the distorted results, especially near local inhomogeneities (e.g. Ley-Cooper et al.2015; Persova et al.2015b, Ullmann et al.2016). In induction logging, this problem is important, e.g. when investigating the horizontal wells, and if the geometry and electrophysical properties of the borehole environment are locally and drastically changed or it is necessary to correctly define the intersection of the fault by the borehole. In these cases, for adequately recovering the medium, it is necessary to use the multidimensional approaches based on accurate (without significant simplifications) mathematical models of the electromagnetic field. Therefore, the increasingly important problem is the development of such methods of electromagnetic fields calculation that, on the one hand, would make it possible to model the field in a complex medium with the required adequacy (with the use of mathematical models without simplifications) and, on the other hand, would be computationally low-cost and convenient to use. For multidimensional electromagnetic fields modelling in complex 3-D media, among the mesh-based numerical methods the following are more frequently used: the finite-difference (FD) method (Grayver et al.2013; Jaysaval et al.2015; Shantsev & Maaø 2015, etc.), the finite volume method (Haber & Ruthotto 2014; Haber & Schwarzbach 2014; Jahandari & Farquharson 2014; Jahandari & Farquharson 2015, etc.) and the finite-element (FE) method (Farquharson & Miensopust 2011; Ansari & Farquharson 2014; Börner et al.2015; Persova et al.2015a; Kordy et al.2016a, etc.). Besides, for modelling the multidimensional electromagnetic fields including the solution of 3-D induction logging problems, the integral equation method is applied quite successfully (Fang et al.2003; Nie et al.2013; Yang et al.2017, etc.). Its computational effectiveness is especially high if the volume of 3-D inhomogeneities is small. When solving the multidimensional electromagnetic problems with the use of mesh-based numerical methods, a significant lowering of computational costs can be achieved by means of the following approaches. First, this is the primary-secondary field approach which uses the separate calculation of the horizontally layered medium field (that is called the primary field). Among the first papers, which considered this approach for FD and FE methods, one should mention the papers (Newman & Alumbaugh 1995; Soloveichik et al.1997; Soloveichik et al.1998; Badea et al.2001). Nowadays, this approach continues to be actively applied (e.g. Mukherjee & Everett 2011; Zaslavsky et al.2011; Da Silva et al.2012; Grayver et al.2013). The more similar the investigated medium to the horizontally layered one, which can be taken different for different logging tool positions, the higher the effectiveness of the primary-secondary field approach. In this case, for achieving the required accuracy of the summary (total) field calculation, when the anomalous (secondary) field is calculated, the number of nodes in the 3-D mesh can be decreased by several times, and the computational time of one 3-D problem solution can be decreased by an order or more. Second, a significant lowering of the computational costs can be achieved if we use the unstructured meshes (e.g. Lelièvre & Farquharson 2013; Wang et al.2013; Um et al.2015) including the non-conforming ones with local refinements (e.g. Grayver & Bürg 2014; Persova et al.2014b; Grayver 2015). In this case, in comparison with the regular meshes, the number of elements in the mesh can be decreased by several times (without the loss of the accuracy of the 3-D problem solution), and this allows the additional reduction of the computational costs by an order. Third, if the number of the logging tool positions is big, a significant lowering of the computational costs can be achieved if we use the direct solvers for the finite-difference or FE equations system (Da Silva et al.2012; Oldenburg et al.2013; Chung et al.2014; Haber & Schwarzbach 2014; Puzyrev et al.2016) (e.g. PARDISO solver (Schenk & G̈artner 2004)). In this case, the matrix factorization, which is the most computational cost procedure in calculating a 3-D field, is performed only once for several logging tool positions simultaneously, and the computational time of 3-D problems for all logging tool positions joined in one relatively small group is close to the computational time of the 3-D problem for one tool position. However, it is necessary to take into account that the computational time for one tool position when the direct solver is used can be much higher in comparison with the time when the iteration solver is used. In this paper, we consider the approach to calculating the 3-D field in the induction logging problems, in which all of the three above methods are used for lowering the computational costs. Additionally, we propose: a new method of FE calculation of the primary field in a horizontally layered medium, which is effective when it is necessary to calculate the primary field at many points simultaneously (that is necessary for calculating the secondary field); the method of 3-D electromagnetic field FE calculation on irregular non-conforming meshes with the construction of conforming basis functions from the non-conforming ones and corresponding procedure of assembling matrices of FEs; the method of grouping the logging tool positions (in order to use the direct solver when calculating the 3-D field for all positions in one group) with the automatic choice of the 1-D background medium for calculating the primary fields and automatic generation of the irregular 3-D mesh with hanging nodes for one group of logging tool positions. 2 MATHEMATICAL MODELS 2.1 The model describing the electromagnetic field in a 3-D anisotropic medium We represent a total electric field in the form $${{{\bf \vec{E}}}^t} = {{{\bf \vec{E}}}^p} + {{{\bf \vec{E}}}^s}$$ where $${{{\bf \vec{E}}}^p}$$ is the primary electric field induced by a source in a horizontally layered (1-D background) medium, and $${{{\bf \vec{E}}}^s}$$ is the 3-D secondary electric field determined by the impact of 3-D inhomogeneities such as layers pinchout, faults, non-parallelism of layers, etc. We can seek $${{{\bf \vec{E}}}^s}$$ by solving the following vector differential equation:   \begin{equation} \frac{1}{{{\mu _0}}}\nabla \times \nabla \times {{{\bf \vec{E}}}^s} + i\omega \sigma {{{\bf \vec{E}}}^s} - {\omega ^2}\varepsilon {{{\bf \vec{E}}}^s} = - i\omega \left( {\sigma - {\sigma ^p}} \right){{{\bf \vec{E}}}^p} + {\omega ^2}\varepsilon {{{\bf \vec{E}}}^p}, \end{equation} (1)where μ0 denotes the vacuum magnetic permeability, σ denotes the 3 × 3 conductivity tensor in a 3-D medium, σp denotes the 3 × 3 conductivity tensor in a horizontally layered (1-D background) medium (for which $${{{\bf \vec{E}}}^p}$$ is calculated), ε(x, y, z) denotes the dielectric permeability of the 3-D medium, ω = 2πν denotes the angular frequency (ν is the electromagnetic field frequency), and $$i = \sqrt { - 1} $$. The components σij(x, y, z) of the tensor σ and components $$\sigma _{ii}^p( z )$$ of the tensor σp are the piecewise constant functions of coordinates. The components σii are not equal to $$\sigma _{ii}^p$$ only in the subdomains where 3-D inhomogeneities are located. When i ≠ j, σij are not equal to zero only in the subdomains where 3-D anisotropic inhomogeneities are located, and the anisotropic axes of these inhomogeneities do not coincide with the coordinate axes. Dielectric permeability ε can also be defined as a tensor. If thin layers in the anisotropic medium are not parallel to Cartesian planes (i.e. anisotropic axes do not coincide with coordinate axes), the conductivity tensor can be obtained by the relation   \begin{equation} \sigma = {{{\bf M}}^{{\bf T}}}\tilde{\sigma }{{\bf M}}, \end{equation} (2)where M is the rotation matrix, $$\tilde{\sigma }$$ is the diagonal tensor with the components corresponding to the lateral and vertical conductivities (in the local coordinate system related with anisotropic axes), the superscript ‘T’ denotes transposition. In a similar way, the tensor of dielectric permeability can be defined in the anisotropic medium. Eq. (1) defines the secondary field for the case when the primary field is calculated without the account of displacement currents. In general, in the induction logging problems, when electromagnetic field is calculated for the frequency ν < 200 kHz, the displacement currents effects are minor, and the terms with ε can be excluded from eq. (1) (i.e. ε is assumed to be zero). If the frequency ν > 200 kHz, the displacement currents effect on the calculated electromagnetic field can be significant. In this case, the displacement currents are to be taken into account in calculating the secondary field. In induction logging measurements, the transmitters (induction coils) with different orientation can be used. Most often the coils are coaxial with the borehole. But sometimes the tools with several (two or three) orthogonal coils (the axis of one of these coils coincides with the tool axis) are considered. For computing $${{{\bf \vec{E}}}^p}$$ induced by an arbitrary oriented coil in the horizontally layered background medium, it is required to compute two fields for two coil orientations: in the first case the coil axis is perpendicular to the layers, in the second one the coil axis is parallel to the layers. In general, for calculating $${{{\bf \vec{E}}}^p}$$, analytical or semi-analytical methods can be used. The first and the simplest is the approach where the field of homogeneous space is taken as $${{{\bf \vec{E}}}^p}$$. This field can be calculated by the analytical method. The disadvantage of this approach is the following. Even if the conductivity of the homogeneous space is taken equal to the conductivity of the layer in which the transmitter is located, it is impossible to reach a comparable effectiveness that can be achieved if we take the horizontally layered medium field as $${{{\bf \vec{E}}}^p}$$. This is because, in the first case, the layers (as 3-D inhomogeneities) are the sources of the secondary field. Therefore, for reaching the required accuracy, it is necessary to use sufficiently fine meshes (that will lead to much higher computational costs). Until recently, semi-analytical methods have been developed quite intensively (Key 2009). They are very effective if there is a need to compute the field at a relatively moderate number of points (corresponding to receiver locations). But if it is required to compute the field values at a great number of points, the computational efficiency of these methods is not so high. For implementing the method of calculating the 3-D electromagnetic field when the horizontally layered medium is taken as a background medium, it is required to compute the primary field $${{{\bf \vec{E}}}^p}$$ at a great number of points (cells) in subdomains (of the computational domain) where 3-D inhomogeneities are located. Even if the optimized irregular non-conforming meshes are used, the primary field needs to be calculated at tens of thousands of points. Consequently, in this case the approaches which allow us to compute the horizontally layered medium field in all the investigated volume simultaneously are more advantageous (Persova et al.2011; Mogilatov et al.2016). Therefore, for calculating the primary field $${{{\bf \vec{E}}}^p}$$, we propose another approach based on the use of the FE method. 2.2 Calculating $${{{\bf \vec{E}}}^p}$$ in a horizontally layered background medium We represent $${{{\bf \vec{E}}}^p}$$ as a sum of solutions of several axisymmetrical problems. When describing mathematical models in the cylindrical coordinate system, we use the diagonal 2 × 2 tensor σp, rz(z) with the components $$\sigma _r^{p,rz} = \sigma _{11}^p = \sigma _{22}^p$$ and $$\sigma _z^{p,rz} = \sigma _{33}^p$$ (where $$\sigma _{ii}^p$$ are the diagonal components of the tensor σp in eq. 1). At first, we consider the mathematical model for computing the electrical field $${{{\bf \vec{E}}}^{p,1}}$$ which is induced by a ‘horizontal’ coil. This coil lies in plane which is parallel to the layers of a horizontally layered background medium (and, obviously, the coil axis is perpendicular to the layers borders). $${{{\bf \vec{E}}}^{p,1}}$$ induced by this coil is axisymmetrical and has only one non-zero component $$E_\varphi ^{}( {r,z} )$$ in cylindrical coordinates {r, φ, z}. This component ($$E_\varphi ^{}( {r,z} )$$) can be found by solving the equation   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta E_\varphi ^{} + \frac{1}{{{\mu _0}{r^2}}}E_\varphi ^{} + i\omega \sigma _r^{p,rz}E_\varphi ^{} = - i\omega I{\delta ^{{J_\varphi }}}, \end{equation} (3)where I is the value of the coil current, and $${\delta ^{{J_\varphi }}}$$ is the delta-function concentrated on the circle r = rδ, z = zδ corresponding to the coil contour ($$I{\delta ^{{J_\varphi }}}$$ is the source current density Jφ). The components of the vectors $${{{\bf \vec{E}}}^{p,1}}$$ and $${{{\bf \vec{B}}}^{p,1}}$$ can be obtained using $$E_\varphi ^{}$$ by the following relations:   \begin{equation} E_x^{p,1}\left( {x,y,z} \right) = - E_\varphi ^{}\left( {r,z} \right)\frac{{y - {y^c}}}{r},\quad E_y^{p,1}\left( {x,y,z} \right) = E_\varphi ^{}\left( {r,z} \right)\frac{{x - {x^c}}}{r},\quad E_z^p = 0, \end{equation} (4)  \begin{equation*} B_x^{p,1}\left( {x,y,z} \right) = \frac{1}{{i\omega }}\frac{{\partial {E_\varphi }\left( {r,z} \right)}}{{\partial z}}\frac{{x - {x^c}}}{r},\quad B_y^{p,1}\left( {x,y,z} \right) = \frac{1}{{i\omega }}\frac{{\partial {E_\varphi }\left( {r,z} \right)}}{{\partial z}}\frac{{y - {y^c}}}{r}, \end{equation*}   \begin{equation} B_z^{p,1}\left( {x,y,z} \right) = - \frac{1}{{i\omega }}\left( {\frac{{{E_\varphi }\left( {r,z} \right)}}{r} + \frac{{\partial {E_\varphi }\left( {r,z} \right)}}{{\partial r}}} \right), \end{equation} (5)where xc and yc are the coordinates of the coil centre, and $$r = \sqrt {{{( {x - {x^c}} )}^2} + {{( {y - {y^c}} )}^2}} $$. The electrical field $${{{\bf \vec{E}}}^{p,2}}$$ induced by a ‘vertical’ coil (its axis is parallel to the borders of the layers of the horizontally layered medium) is 3-D even for a 1-D (horizontally layered) medium, that is, the components of $${{{\bf \vec{E}}}^{p,2}}$$ depend on all three space coordinates. In order to compute this field, we apply the approach similar to the one considered by Persova et al. (2011) for modelling the field induced by the horizontal electrical dipole. It is based on the representation of the field, which is in fact a 3-D field in the 1-D medium, as a sum of the fields each of which can be obtained by the solution of the problem of a less dimension. However, when we use the approach considered by Persova et al. (2011) there appears a problem (this is particularly true for harmonic modes) which is related to the fact that the total field is formed by the components with the opposite signs the moduli of which are close to each other. This makes it necessary to increase significantly the accuracy of the components calculation that, in its turn, leads to the significant increase in computational costs. Therefore, in this paper we propose the modification of the approach described by Persova et al. (2011) that makes it possible to obtain $${{{\bf \vec{E}}}^{p,2}}$$ (the field induced by the ‘vertical’ coil) as a sum of several axisymmetrical fields without mutually subtracting components (Soloveichik et al.2014; Persova et al.2014a). Each of these fields (forming $${{{\bf \vec{E}}}^{p,2}}$$) is induced by a special source and can be obtained by solving the corresponding axisymmetrical problem. At first we note that, in the problem considered, the sizes of the transmitter coil are small relative to the distance to the receivers. Therefore, we can take any geometrical form of coil. We represent the ‘vertical’ coil as a rectangular (square) frame. In this case the total source (frame) can be represented as a sum of sources shown in Fig. 1. This representation allows us to obtain the electromagnetic field in a horizontally layered background medium by means of solutions of several axisymmetrical problems. Figure 1. View largeDownload slide The view of (a) non-axisymmetric source ‘I1’ and (b) two axisymmetrical sources ‘I2’; (c) the current of the upper part of the source ‘I1’ in plane XY; (d) the currents of one of NGR sources in plane XY; (e) currents configuration which is obtained when the negative source ‘I2’ is added to the source ‘I1’; (f) currents configuration which is obtained when the positive source ‘I2’ is added to the source shown in panel (e) (as a result, the ‘vertical’ rectangular current frame is obtained); (g) the approximation of the horizontal ‘current lines’ by the point sources located in the centres of subintervals of these ‘current lines’ (the length of subintervals is $$h_L^{}$$). Dashed light-green arrows show the direction of the surface radial currents flowing in plane XY. Solid black lines show the currents flowing in the ‘vertical’ current frame. Red points show the positions of the point sources approximating the horizontal ‘current lines’. Figure 1. View largeDownload slide The view of (a) non-axisymmetric source ‘I1’ and (b) two axisymmetrical sources ‘I2’; (c) the current of the upper part of the source ‘I1’ in plane XY; (d) the currents of one of NGR sources in plane XY; (e) currents configuration which is obtained when the negative source ‘I2’ is added to the source ‘I1’; (f) currents configuration which is obtained when the positive source ‘I2’ is added to the source shown in panel (e) (as a result, the ‘vertical’ rectangular current frame is obtained); (g) the approximation of the horizontal ‘current lines’ by the point sources located in the centres of subintervals of these ‘current lines’ (the length of subintervals is $$h_L^{}$$). Dashed light-green arrows show the direction of the surface radial currents flowing in plane XY. Solid black lines show the currents flowing in the ‘vertical’ current frame. Red points show the positions of the point sources approximating the horizontal ‘current lines’. The rectangular ‘vertical’ current frame is represented as a sum of source ‘I1’ shown in Fig. 1(a) (we consider the top and bottom currents as the parts of one source ‘I1’) and two sources ‘I2’ shown in Fig. 1(b). The source ‘I1’ consists of four surface radial currents, which flow in plane XY (in Fig. 1, their direction is shown by dashed light-green arrows), and two horizontal ‘current lines’ (in Fig. 1, the currents in them are shown by solid black lines). We consider the left and right currents as different sources of the type ‘I2’, the former is negative and the latter is positive with regard to the direction of the axes shown in Fig. 1 (in this case we direct the z-axis down because in all examples presented in this paper, the direction of the coordinate axes is the same). Each of the sources ‘I2’ consists of the vertical ‘current line’ (the current in which shown by solid black line) and two surface radial currents (their direction is shown by dashed light-green arrows). Joining the top and bottom currents in one source ‘I1’ and separating the left and right currents into two sources of the type ‘I2’ are convenient for further reduction of the fields induced by these sources to the corresponding axisymmetrical fields. The Z-coordinate of the upper part of source ‘I1’ coincides with the Z-coordinate of the upper frame border, and the Z-coordinate of the lower part of the source ‘I1’ (with the opposite sign of the current) coincides with the Z-coordinate of the lower frame border. The axis of one of the sources ‘I2’ coincides with the left frame border, and the axis of another source ‘I2’ (with an opposite sign of the current) coincides with the right frame border. Convergent and divergent surface radial currents are the parts of these sources. They provide zero current divergence in sources ‘I1’ and ‘I2’. Further, we will refer to as ‘non-grounded radial sources’ (NGR sources) for each of these surface radial currents converging into the point Pi or diverging from it. Fig. 1(c) shows the current of the upper part of the source ‘I1’ in plane XY, and Fig. 1(d) shows the currents of one of NGR sources in plane XY. For each of NGR sources, the same current passes through any circle C the centre of which coincides with the NGR source centre (see Fig. 1d, the integrals along circles С1 and С2 of surface radial currents of the NGR source are equal). Figs 1(e) and (f) demonstrate how the current in the ‘vertical’ frame is composed from the source ‘I1’ and two sources ‘I2’. In Fig. 1(e) the negative source ‘I2’ with the axis, which coincides with the left side of the ‘vertical’ frame, is added to the source ‘I1’. The radial currents of the source ‘I2’ compensate the radial currents of the source ‘I1’ in the left top and left bottom angles (these currents are equal and are in the opposite directions). When the positive source ‘I2’ (with the axis coinciding with the right side of the ‘vertical’ frame) is added to the source shown in Fig. 1(e), the radial currents in the right top and right bottom angles are compensated. As a result, we obtain the source with the current in the ‘vertical’ frame shown in Fig. 1(f). We would like to reiterate that each of the sources ‘I2’ generates the axisymmetrical field (in the coordinate system in which the z-axis coincides with the vertical current of the source). The mathematical model for calculating the axisymmetrical field generated by the source ‘I2’ is presented below (see eq. 21–23). We first consider how the field induced by the non-axisymmetric source ‘I1’ can be calculated by solving several axisymmetrical problems. We introduce the vector potential $${{\bf \vec{A}}}$$ that is related with the magnetic induction and electric field by the relationships $${{\bf \vec{B}}} = {\rm{rot}}{{\bf \vec{A}}}$$ and $${{\bf \vec{E}}} = - \frac{{\partial {{\bf \vec{A}}}}}{{\partial t}} = - i\omega {{\bf \vec{A}}}$$. Then, the vector potential $${{{\bf \vec{A}}}^{{\rm{I1}}}}$$ defining the electromagnetic field induced by source ‘I1’ in a horizontally layered medium (with conductivity σp) can be found by solving the equation   \begin{equation} \frac{1}{{{\mu _0}}}\nabla \times \nabla \times {{{\bf \vec{A}}}^{{\rm{I1}}}} + i\omega {\sigma ^p}{{{\bf \vec{A}}}^{{\rm{I1}}}} = {{{\bf \vec{J}}}^{{\rm{I1}}}}, \end{equation} (6)where $${{{\bf \vec{J}}}^{{\rm{I1}}}}$$ is the density of the currents defining source ‘I1’. We define this source in such a way that, first, its currents have only horizontal non-zero components, and, second, for any bounded subdomain Ωf, the value of the current entering Ωf is equal to the value of the outward current. Therefore, for the current density $${{{\bf \vec{J}}}^{{\rm{I1}}}}$$, the ratio   \begin{equation} \!\nabla \cdot {{{\bf \vec{J}}}^{{\rm{I1}}}} = 0 \end{equation} (7)is true and $$J_z^{{\rm{I1}}} \equiv 0$$. It would now be easy to show that eq. (6) is satisfied for the vector potential $${{{\bf \vec{A}}}^{{\rm{I1}}}}$$ whose z-component is equal to zero and the components $$A_x^{{\rm{I1}}}$$ and $$A_y^{{\rm{I1}}}$$ satisfy the equations   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta A_x^{{\rm{I1}}} + i\omega \sigma _{11}^pA_x^{{\rm{I1}}} = J_x^{{\rm{I1}}}, \end{equation} (8)  \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta A_y^{{\rm{I1}}} + i\omega \sigma _{22}^pA_y^{{\rm{I1}}} = J_y^{{\rm{I1}}}, \end{equation} (9)where Δ denotes the Laplacian, and $$\sigma _{11}^p = \sigma _{22}^p = \sigma _r^{p,rz}(z)$$. We will prove it. For this purpose, we differentiate eq. (8) with respect to x, and eq. (9) with respect to y (taking into account that $$\sigma _{11}^p$$ and $$\sigma _{22}^p$$ do not depend on x and y and equal to $$\sigma _r^{p,rz}$$). Changing the order of differentiation in the first summands of the left-hand sides of these equations and summing them, we obtain   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta \left( {\frac{{\partial\! A_x^{{\rm{I1}}}}}{{\partial x}} + \frac{{\partial\! A_y^{{\rm{I1}}}}}{{\partial y}}} \right) + i\omega \sigma _r^{p,rz}\left( {\frac{{\partial\! A_x^{{\rm{I1}}}}}{{\partial x}} + \frac{{\partial\! A_y^{{\rm{I1}}}}}{{\partial y}}} \right) = \left( {\frac{{\partial J_x^{{\rm{I1}}}}}{{\partial x}} + \frac{{\partial J_y^{{\rm{I1}}}}}{{\partial y}}} \right). \end{equation} (10)Taking into account that $$A_z^{{\rm{I1}}} \equiv 0$$ and $$J_z^{{\rm{I1}}} \equiv 0$$ (and, consequently, $$\frac{{\partial\! A_z^{{\rm{I1}}}}}{{\partial z}} = 0$$ and $$\frac{{\partial J_z^{{\rm{I1}}}}}{{\partial z}} = 0$$), we obtain   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta ( {\nabla \cdot {{{{\bf \vec{A}}}}^{{\rm{I1}}}}}) + i\omega \sigma _r^{p,rz}\left( {\nabla \cdot {{{{\bf \vec{A}}}}^{{\rm{I1}}}}} \right) = \nabla \cdot {{{\bf \vec{J}}}^{{\rm{I1}}}}. \end{equation} (11)As far as eq. (7) is true for $${{{\bf \vec{J}}}^{{\rm{I1}}}}$$, the function $$u = \nabla \cdot {{{\bf \vec{A}}}^{{\rm{I1}}}}$$ satisfies the homogeneous equation   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta u + i\omega \sigma _r^{p,rz}u = 0 \end{equation} (12)for which the function u = 0 is the solution. Thus, the divergence of vector potential $${{{\bf \vec{A}}}^{{\rm{I1}}}}$$ (whose x- and y-components satisfy eqs (8) and (9), and $$A_z^{{\rm{I1}}} \equiv 0$$) is equal to zero. Taking into account $$A_z^{{\rm{I1}}} = 0$$, the equation system (8), (9) can be written in the vector form   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta {{{\bf \vec{A}}}^{{\rm{I1}}}} + i\omega {\sigma ^p}{{{\bf \vec{A}}}^{{\rm{I1}}}} = {{{\bf \vec{J}}}^{{\rm{I1}}}}, \end{equation} (13)where Δ denotes the vector Laplacian. As far as $$\nabla \cdot {{{\bf \vec{A}}}^{{\rm{I1}}}} = 0$$ (we have just shown it), and taking into account the identical relation $$\nabla \times \nabla \times {{{\bf \vec{A}}}^{{\rm{I1}}}} \equiv - \Delta {{{\bf \vec{A}}}^{{\rm{I1}}}} + \nabla \nabla \cdot {{{\bf \vec{A}}}^{{\rm{I1}}}}$$, the vector potential $${{{\bf \vec{A}}}^{{\rm{I1}}}}$$, which is obtained by solving eqs (8) and (9) (or by solving vector eq. 13), satisfies eq. (6). Now we show how the solution of eq. (13) could be obtained by solving several axisymmetrical problems. For this purpose, we divide the source ‘I1’ into the following parts: four NGR sources and two ‘current lines’. By ‘current line’ we mean the source which is concentrated on one of the horizontal borders P1P2 or P3P4 of the ‘vertical’ current frame (see Fig. 1). The density of such sources can be defined via delta-function δL, which is concentrated on the line and defined in the following way: the integral of δL taken over an arbitrary volume Ωf is equal to the length of the L line segment, which lies inside the volume Ωf. Thus, we represent right-hand side of vector eq. (13) as a sum of three summands:   \begin{equation} {{{\bf \vec{J}}}^{{\rm{I1}}}} = {{{\bf \vec{J}}}^{{\rm{I1,NGR1}}}} + {{{\bf \vec{J}}}^{{\rm{I1,NGR2}}}} + {{{\bf \vec{J}}}^{{\rm{I1,hor}}}}, \end{equation} (14)where $${{{\bf \vec{J}}}^{{\rm{I1,NGR1}}}}$$ is the density of surface currents in two (axial) NGR sources, while the currents of one of these sources converge into the point P1, and currents of another source diverge from the point P3 (see Fig. 1a); $${{{\bf \vec{J}}}^{{\rm{I1,NGR2}}}}$$ is the density of surface currents in other two (axial) NGR sources (currents of one of these sources diverge from the point P2, and currents of another source converge to the point P4); $${{{\bf \vec{J}}}^{{\rm{I1,hor}}}} = {{{\bf \vec{J}}}^{{\rm{I1,}}{{\rm{P}}_1}{{\rm{P}}_2}}} - {{{\bf \vec{J}}}^{{\rm{I1,}}{{\rm{P}}_3}{{\rm{P}}_4}}}$$ is the density of currents in the horizontal ‘current lines’ P1P2 and P3P4. Then, the solution of eq. (13) can be found as a sum of solutions of three equations (they are in fact independent)   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta {{{\bf \vec{A}}}^{{\rm{I1,NGR1}}}} + i\omega {\sigma ^p}{{{\bf \vec{A}}}^{{\rm{I1,NGR1}}}} = {{{\bf \vec{J}}}^{{\rm{I1,NGR1}}}}, \end{equation} (15)  \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta {{{\bf \vec{A}}}^{{\rm{I1,NGR2}}}} + i\omega {\sigma ^p}{{{\bf \vec{A}}}^{{\rm{I1,NGR2}}}} = {{{\bf \vec{J}}}^{{\rm{I1,NGR2}}}}, \end{equation} (16)  \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta {{{\bf \vec{A}}}^{{\rm{I1,hor}}}} + i\omega {\sigma ^p}{{{\bf \vec{A}}}^{{\rm{I1,hor}}}} = {{{\bf \vec{J}}}^{{\rm{I1,}}{{\rm{P}}_1}{{\rm{P}}_2}}} - {{{\bf \vec{J}}}^{{\rm{I1,}}{{\rm{P}}_3}{{\rm{P}}_4}}}. \end{equation} (17)Right-hand sides of eqs (15) and (16) are the functions with axial symmetry (the z-axis passes through the points P1,P3 or P2,P4, respectively). In the corresponding cylindrical coordinate system, the functions $${{{\bf \vec{J}}}^{{\rm{I1,NGR1}}}}$$ and $${{{\bf \vec{J}}}^{{\rm{I1,NGR2}}}}$$ of the right-hand sides of these equations have the only non-zero radial component: $$J_r^{{\rm{I1,NGR1}}} = \frac{I}{{2\pi r}}( {{\delta ^{{{\hat{\Gamma }}_2}}} - {\delta ^{{{\hat{\Gamma }}_1}}}} )$$ and $$J_r^{{\rm{I1,NGR2}}} = \frac{I}{{2\pi r}}( {{\delta ^{{{\hat{\Gamma }}_1}}} - {\delta ^{{{\hat{\Gamma }}_2}}}} )$$, where I is the value of current in the ‘vertical’ frame, $${\delta ^{{{\hat{\Gamma }}_1}}}$$ and $${\delta ^{{{\hat{\Gamma }}_2}}}$$ are the delta-functions concentrated on the surfaces $${\hat{\Gamma }_1}$$ and $${\hat{\Gamma }_2}$$ (i.e. on the planes z = z1 and z = z2 corresponding to top and bottom frame borders, see Fig. 1) where the radial currents of the top and bottom NGR sources flow. Thus, the solution of both eqs (15) and (16) can be found as a result of solving one axisymmetrical problem: find the solution $$A_r^{}( {r,z} )$$ of equation   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta A_r^{} + \frac{1}{{{\mu _0}{r^2}}}A_r^{} + i\omega \sigma _r^{p,rz}A_r^{} = \frac{I}{{2\pi r}}( {{\delta ^{{{\hat{\Gamma }}_1}}} - {\delta ^{{{\hat{\Gamma }}_2}}}} ), \end{equation} (18)where $$\sigma _r^{p,rz}$$ is the lateral conductivity of a horizontally layered medium. The field $${{{\bf \vec{A}}}^{{\rm{I1,hor}}}}$$ defined by eq. (17) is not axisymmetrical. The vector function $${{{\bf \vec{J}}}^{{\rm{I1,}}{{\rm{P}}_1}{{\rm{P}}_2}}} - {{{\bf \vec{J}}}^{{\rm{I1,}}{{\rm{P}}_3}{{\rm{P}}_4}}}$$ defining the right-hand side of this equation has only one non-zero component depending on the orientation of the ‘vertical’ current frame: if the frame is oriented in the way shown in Fig. 1, then the y-component is the only non-zero component of the right-hand side vector of eq. (17) and, consequently, the vector function $${{{\bf \vec{A}}}^{{\rm{I1,hor}}}}$$ has the only non-zero component $$A_y^{{\rm{hor}}}( {x,y,z} )$$. If the ‘vertical’ current frame is oriented so that the lines P1P2 and P3P4 are parallel to the x-axis, then the x-component is the only non-zero component of the right-hand side vector of eq. (17) and, therefore, the vector function $${{{\bf \vec{A}}}^{{\rm{I1,hor}}}}$$ has the only non-zero component $$A_x^{{\rm{hor}}}( {x,y,z} )$$. These components can be found by solving the 3-D scalar equation   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta A_\xi ^{{\rm{hor}}} + i\omega \sigma _r^{p,rz}A_\xi ^{{\rm{hor}}} = I\left( {{\delta ^{{{\rm{P}}_1}{{\rm{P}}_2}}} - {\delta ^{{{\rm{P}}_3}{{\rm{P}}_4}}}} \right), \end{equation} (19)where $$A_\xi ^{{\rm{hor}}}( {x,y,z} )$$ is the vector potential component corresponding to the horizontal line direction (x-direction or y-direction), $${\delta ^{{{\rm{P}}_1}{{\rm{P}}_2}}}$$ and $${\delta ^{{{\rm{P}}_3}{{\rm{P}}_4}}}$$ are delta-functions concentrated on the lines P1P2 and P3P4. The right-hand side of eq. (19) can be approximated with any required accuracy by the set of point sources positioned along the lines P1P2 and P3P4 (see Fig. 1g) and, therefore, the 3-D component $$A_\xi ^{{\rm{hor}}}( {x,y,z} )$$ of the field induced by the sources concentrated on the lines P1P2 and P3P4 can be computed as a sum of axisymmetrical fields $$A_\xi ^{}( {r,z} )$$ induced by the point sources positioned along these lines. The field induced by two point sources (lying one below (above) another) can be calculated by the solution of the equation   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta A_\xi ^{} + i\omega \sigma _r^{p,rz}A_\xi ^{} = I{h_L}\left( {{\delta ^{{z_1}}} - {\delta ^{{z_2}}}} \right), \end{equation} (20)where $$A_\xi ^{}( {r,z} )$$ is the vector potential component corresponding to the ‘current line’ direction (i.e. x- or y-component), $${\delta ^{{z_1}}}$$ and $${\delta ^{{z_2}}}$$ are the delta-functions concentrated in the points (centres of the ‘current line’ subintervals, in which the point sources are located), hL is the length of the ‘current line’ subinterval (see Fig. 1g), and, I as before, is the value of the current in the frame. By summing axisymmetrical fields $$A_\xi ^{}( {r,z} )$$ taking into account the translation of the cylindrical coordinate system centres into the corresponding points of one of the ‘current lines’ of source ‘I1’, we obtain the 3-D distribution of one component of the vector potential, which describes the field of the ‘current line’. Thus, we can calculate the field induced by the non-axisymmetric source ‘I1’ by summing axisymmetrical fields taking into account the translation of the cylindrical coordinate system centres into corresponding points. The formulae for calculating the electric and magnetic field will be presented below (we take into account the fact that, when summing the components of the fields induced by separate parts of the source, the radial currents are compensated, and, therefore, the fields induced by these radial currents are not added to the resulting electric and magnetic field). In order to compute the axisymmetrical field induced by the source of type ‘I2’, we use the following mathematical model. Taking into account that the divergence of current density $${{{\bf \vec{J}}}^{{\rm{I2}}}} = {( {{J_r},{J_z}} )^{\rm{T}}}$$ is equal to zero, the electromagnetic field induced by each source ‘I2’ can be obtained by the solution of the differential equation system   \begin{equation} \! - \frac{{\rm{1}}}{{{\mu _0}}}\Delta A_r^{} + \frac{{\rm{1}}}{{{\mu _0}{r^2}}}A_r^{} + \sigma _r^{p,rz}\left( {i\omega A_r^{} + \frac{{\partial V}}{{\partial r}}} \right) = {J_r}, \end{equation} (21)  \begin{equation} \! - \frac{{\rm{1}}}{{{\mu _0}}}\Delta A_z^{} + \sigma _z^{p,rz}\left( {i\omega A_z^{} + \frac{{\partial V}}{{\partial z}}} \right) = {J_z}, \end{equation} (22)  \begin{equation} \! - \nabla \cdot \left( {{\sigma ^{p,rz}}\nabla V} \right) - i\omega \nabla \cdot \left( {{\sigma ^{p,rz}}{{{{\bf \vec{A}}}}^{{\rm{I2}}}}} \right) = 0, \end{equation} (23)where $${{{\bf \vec{A}}}^{{\rm{I2}}}} = {( {{A_r}( {r,z} ),{A_z}( {r,z} )} )^{\rm{T}}}$$ is the vector potential in the cylindrical coordinate system, V = V(r, z) is the scalar potential, and Δ = ∇ · ∇ is the scalar Laplacian. We represent the r-component of the vector potential $${{{\bf \vec{A}}}^{{\rm{I2}}}}$$ as $$A_r^{} = A_r^* + \tilde{A}_r^{}$$ where $$A_r^*$$ is the solution of equation in the form of (18). Then system (21)–(23) takes the form   \begin{equation} \! - \frac{{\rm{1}}}{{{\mu _0}}}\Delta \tilde{A}_r^{} + \frac{{\rm{1}}}{{{\mu _0}{r^2}}}\tilde{A}_r^{} + \sigma _r^{p,rz}\left( {i\omega \tilde{A}_r^{} + \frac{{\partial V}}{{\partial r}}} \right) = 0, \end{equation} (24)  \begin{equation} \! - \frac{{\rm{1}}}{{{\mu _0}}}\Delta A_z^{} + \sigma _z^{p,rz}\left( {i\omega A_z^{} + \frac{{\partial V}}{{\partial z}}} \right) = {J_z}, \end{equation} (25)  \begin{equation} \! - \nabla \cdot \left( {{\sigma ^{p,rz}}\nabla V} \right) - i\omega \nabla \cdot \left( {{\sigma ^{p,rz}}{{( {\tilde{A}_r^{},A_z^{}} )}^{\rm{T}}}} \right) = i\omega \nabla \cdot \left( {{\sigma ^{p,rz}}{{( {A_r^*,0} )}^{\rm{T}}}} \right). \end{equation} (26)The field $$A_r^*$$ in the right-hand side of eq. (26) is fully identical with $$A_r^{}$$ calculated for the source ‘I1’ (with the opposite sign) that allows us to solve eq. (18) only once and (what is more important!) to exclude these two mutually subtracting fields when we form the total field of the ‘vertical’ frame. Thus, in order to compute the field of the ‘vertical’ frame, it is necessary to solve three problems: eq. (18), (20) and equation system (24)–(26). Besides, the proposed scheme can be further simplified (that will increase its computational effectiveness) because there is no need to solve eq. (18). We define function u as $$u = \nabla \cdot {(A_r^*,0)^{\rm{T}}}$$. Differentiating eq. (18) properly we obtain the equation   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta u + i\omega \sigma _r^{p,rz}u = I\left( {{\delta ^{P,1}} - {\delta ^{P,2}}} \right), \end{equation} (27)which is identical to eq. (20). Equation system (24)–(26) takes the form   \begin{equation} \! - \frac{{\rm{1}}}{{{\mu _0}}}\Delta \tilde{A}_r^{} + \frac{{\rm{1}}}{{{\mu _0}{r^2}}}\tilde{A}_r^{} + \sigma _r^{p,rz}\left( {i\omega \tilde{A}_r^{} + \frac{{\partial V}}{{\partial r}}} \right) = 0, \end{equation} (28)  \begin{equation} \! - \frac{{\rm{1}}}{{{\mu _0}}}\Delta A_z^{} + \sigma _z^{p,rz}\left( {i\omega A_z^{} + \frac{{\partial V}}{{\partial z}}} \right) = {J_z}, \end{equation} (29)  \begin{equation} \! - \nabla \cdot \left( {{\sigma ^{p,rz}}\nabla V} \right) - i\omega \nabla \cdot \left( {{\sigma ^{p,rz}}{{( {\tilde{A}_r^{},A_z^{}} )}^{\rm{T}}}} \right) = i\omega \sigma _r^{p,rz}u. \end{equation} (30)Finally, in order to obtain the field of the ‘vertical’ frame, it is necessary to solve two problems: eq. (20) and equation system (28)–(30). Further, we consider how to calculate the components of $${{{\bf \vec{E}}}^{p,2}}$$ using the solution of eq. (20) and equation system (28)–(30). In the case when the ‘vertical’ frame axis is parallel to the x-axis (i.e. the coordinates of its vertexes are (xc, y1, z1), (xc, y2, z1), (xc, y2, z2), and (xc, y1, z2)), the values of the components of $${{{\bf \vec{E}}}^{p,2}}$$ at the point (x, y, z) can be calculated with the help of the following formulae:   \begin{equation} E_x^{p,2}\left( {x,y,z} \right) = \left( {i\omega \tilde{A}_r^{}\left( {{r_1},z} \right) + \frac{{\partial V\left( {{r_1},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {x - {x_c}} \right)}}{{{r_1}}} - \left( {i\omega \tilde{A}_r^{}\left( {{r_2},z} \right) + \frac{{\partial V\left( {{r_2},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {x - {x_c}} \right)}}{{{r_2}}}, \end{equation} (31)  \begin{equation} E_y^{p,2}\left( {x,y,z} \right) = - i\omega A_y^{{\rm{hor}}} + \left( {i\omega \tilde{A}_r^{}\left( {{r_1},z} \right) + \frac{{\partial V\left( {{r_1},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {y - {y_1}} \right)}}{{{r_1}}} - \left( {i\omega \tilde{A}_r^{}\left( {{r_2},z} \right) + \frac{{\partial V\left( {{r_2},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {y - {y_2}} \right)}}{{{r_2}}}, \end{equation} (32)  \begin{equation} E_z^{p,2}\left( {x,y,z} \right) = \left( {i\omega \tilde{A}_z^{}\left( {{r_1},z} \right) + \frac{{\partial V\left( {{r_1},z} \right)}}{{\partial z}}} \right) - \left( {i\omega \tilde{A}_z^{}\left( {{r_2},z} \right) + \frac{{\partial V\left( {{r_2},z} \right)}}{{\partial z}}} \right), \end{equation} (33)where $${r_1} = \sqrt {{{( {x - {x_c}} )}^2} + {{( {y - {y_1}} )}^2}} $$, $${r_2} = \sqrt {{{( {x - {x_c}} )}^2} + {{( {y - {y_2}} )}^2}} $$ and $$A_y^{{\rm{hor}}}( {x,y,z} )$$ is obtained by solving eq. (27) as a result of summing the functions u(r, z) (after recalculating u(r, z) from the cylindrical coordinate system into the Cartesian one), which are calculated for the sources positioned in the centres of subintervals of the horizontal ‘current lines’ of the frame. The components of $${{{\bf \vec{B}}}^{p,2}}$$ at the point (x, y, z) can be calculated with the help of the following formulae:  \begin{equation*} B_x^{p,2}\left( {x,y,z} \right) = - \frac{{\partial\! A_y^{{\rm{hor}}}\left( {r,z} \right)}}{{\partial z}} + \left( {\frac{{\partial {A_r}\left( {{r_1},z} \right)}}{{\partial z}} - \frac{{\partial {A_z}\left( {{r_1},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {y - {y_1}} \right)}}{{{r_1}}} - \left( {\frac{{\partial {A_r}\left( {{r_2},z} \right)}}{{\partial z}} - \frac{{\partial {A_z}\left( {{r_2},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {y - {y_2}} \right)}}{{{r_2}}}, \end{equation*}   \begin{equation*} B_y^{p,2}\left( {x,y,z} \right) = - \left( {\frac{{\partial {A_r}\left( {{r_1},z} \right)}}{{\partial z}} - \frac{{\partial {A_z}\left( {{r_1},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {x - {x^c}} \right)}}{{{r_1}}} + \left( {\frac{{\partial {A_r}\left( {{r_2},z} \right)}}{{\partial z}} - \frac{{\partial {A_z}\left( {{r_2},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {x - {x^c}} \right)}}{{{r_2}}}, \end{equation*}   \begin{equation*} B_z^{p,2}\left( {x,y,z} \right) = \frac{{\partial\! A_y^{{\rm{hor}}}\left( {r,z} \right)}}{{\partial x}}. \end{equation*} If the transmitter coil is oriented so that the angle between its axis and the axis z is equal to α (see Fig. 2a), then $${{{\bf \vec{E}}}^p}$$ can be calculated using $${{{\bf \vec{E}}}^{p,1}}$$ and $${{{\bf \vec{E}}}^{p,2}}$$ by formula  \begin{equation*} {{{\bf \vec{E}}}^p} = {{{\bf \vec{E}}}^{p,1}}\cos \alpha + {{{\bf \vec{E}}}^{p,2}}{\sin \alpha}. \end{equation*} Figure 2. View largeDownload slide (a) The orientation of transmitter whose axis is perpendicular to the y-axis and (b) the rotation of the coordinate system when the transmitter orientation is arbitrary. Figure 2. View largeDownload slide (a) The orientation of transmitter whose axis is perpendicular to the y-axis and (b) the rotation of the coordinate system when the transmitter orientation is arbitrary. For arbitrary transmitter coil orientation, at first, we can calculate $${{{\bf \vec{\tilde{E}}}}^p}$$ by formulae (31)–(33) in the coordinate system $$\{ {\tilde{x},\tilde{y},\tilde{z}} \}$$, in which the coil axis is perpendicular to the $$\tilde{y}$$-axis, and then recalculate $${{{\bf \vec{\tilde{E}}}}^p}$$ into the coordinate system {x, y, z} by formulae $${{{\bf \vec{E}}}^p} = {{\bf M}}{{{\bf \vec{\tilde{E}}}}^p}$$ where M is the rotation matrix defined as follows:   \begin{equation} {{\bf M}} = \left( {\begin{array}{c@{\quad}c@{\quad}c} {\cos \beta }&{ - \sin \beta }&0\\ {\sin \beta }&{\cos \beta }&0\\ 0&0&1 \end{array}} \right). \end{equation} (34)In formula (34), β is the angle of rotation around the axis $$z = \tilde{z}$$ which transfers the coordinate system {x, y} to the coordinate system $$\{ {\tilde{x},\tilde{y}} \}$$(see Fig. 2b). The components of the magnetic field induction are calculated in a similar way. Thus, the proposed approach together with the possibility to take the field of the source in the horizontally layered medium as $${{{\bf \vec{E}}}^p}$$, which is much similar to the realistic 3-D medium in the vicinity of the logging tool location, makes it possible to exclude the appearance of the opposite sign components of the field, the moduli of which are close to each other. 3 FINITE-ELEMENT APPROXIMATION FE solutions of all the equations considered above and describing the secondary field $${{{\bf \vec{E}}}^s}$$ and the components of the primary field (which are axisymmetrical, if the coordinate system is properly translated) are found in the limited 3-D domain Ω (for $${{{\bf \vec{E}}}^s}$$) or 2-D domain Ωrz (for components of $${{{\bf \vec{E}}}^p}$$). On the outer boundary of Ω and Ωrz, we set homogeneous essential boundary conditions. For $${{{\bf \vec{E}}}^s}$$, these conditions represent the equality of $${{{\bf \vec{E}}}^s}$$ components to the zero, which are tangential to the boundary of Ω: $${{{\bf \vec{E}}}^s} \times {{\bf \vec{n}}} = {{\bf \vec{0}}}$$ where $${{\bf \vec{n}}}$$ is the normal to the boundary of Ω (these essential boundary conditions are set when the vector FEM is used). If the field $${{{\bf \vec{E}}}^s}$$ has a symmetry plane, then on this plane, we set either the homogeneous natural boundary condition $$( {\nabla \times {{{{\bf \vec{E}}}}^s}} ) \times {{\bf \vec{n}}} = {{\bf \vec{0}}}$$ (where $${{\bf \vec{n}}}$$ is the normal to the symmetry plane), if the current frame lies in this plane, or the homogeneous essential boundary condition $${{{\bf \vec{E}}}^s} \times {{\bf \vec{n}}} = {{\bf \vec{0}}}$$, if the current frame is sectioned by the symmetry plane. For calculating the axisymmetrical components of the field $${{{\bf \vec{E}}}^p}$$, on the outer boundaries of the computational domain Ωrz, we set homogeneous Dirichlet boundary conditions for all potentials describing the corresponding field. The size of computational domains Ω and Ωrz depends on the frequency and ‘average’ conductivity of the medium. Besides, the size of Ω is much smaller than the size of Ωrz. For example, for the frequency of 200 kHz, Ω = [−20,20] × [−15,15] × [−20,20] m3 (or Ω = [−20,20] × [0,15] × [−20,20] m3 if y = 0 is the symmetry plane of the field), and Ωrz = [0,70] × [−70,70] m2 (if the ‘average’ conductivity of the medium is about 0.05–1 Sm m−1). 3.1 Finite-element approximation for mathematical models used for $${{{\bf \vec{E}}}^p}$$ calculation The equivalent variational formulation for eq. (3), which defines the field of ‘horizontal’ coil (loop), takes the following form:   \begin{equation} \frac{1}{{{\mu _0}}}\left( {\int_{{{\Omega ^{rz}}}}{{\nabla E_\varphi ^{} \cdot \nabla \nu r{\rm{d}}r{\rm{d}}z}} + \int_{{{\Omega ^{rz}}}}{{\frac{1}{r}E_\varphi ^{}\nu {\rm{d}}r{\rm{d}}z}}} \right) + i\omega \int_{{{\Omega ^{rz}}}}{{\sigma _r^{p,rz}E_\varphi ^{}\nu r{\rm{d}}r{\rm{d}}z}} = - i\omega I\nu \left( {{r_\delta },{z_\delta }} \right){r_\delta }, \end{equation} (35)where ν(r, z) is a real test function, Ωrz is the computational domain in the cylindrical coordinate system. The expression in the right-hand side of eq. (35) is obtained as a result of the delta-function $${\delta ^{{J_\varphi }}}$$ (which concentrates on the circle r = rδ, z = zδ) operation on the test function ν(r, z). More specifically, the function $$E_\varphi ^{}$$ satisfying eqs (3) and (35) goes to infinity when approaching the circle r = rδ, z = zδ. The FE solution at the point (rδ, zδ) of the cylindrical domain Ωrz is finite, but together with mesh refinement around this point, the FE solution will grow infinitely. However, at some fixed (non-zero) distance from the point (rδ, zδ), the value of $$E_\varphi ^{}$$ can be calculated with any required accuracy by means of FEM. Therefore, when calculating $${{{\bf \vec{E}}}^{p,1}}$$, this fact must be taken into account, and when forming $${{{\bf \vec{E}}}^p}$$ for the solution of eq. (1), it is not allowed to come to the circle r = rδ, z = zδ closer than some pre-set distance. The solution $$E_\varphi ^{}$$ of eq. (35) is expressed as   \begin{equation} E_\varphi ^{} = \sum\limits_{j = 1}^{{n^{rz}}} {q_j^{{E_\varphi }}{\psi _j}} , \end{equation} (36)where ψj(r, z) are (real) nodal basis functions, nrz is their number (which coincides with the number of nodes in the 2-D mesh). As a result of substituting (36) into (35) and replacing the test function ν with all basis functions ψl in turn, the following system of linear algebraic equations (SLAE) is obtained:   \begin{equation} \!\left( {{{{\bf G}}^{{E_\varphi }}} + i\omega {{{\bf M}}^{{E_\varphi }}}} \right){{{\bf q}}^{{E_\varphi }}} = - i\omega {{{\bf f}}^{{E_\varphi }}}, \end{equation} (37)where the components of the matrices $${{{\bf G}}^{{E_\varphi }}}$$, $${{{\bf M}}^{{E_\varphi }}}$$ and vector $${{{\bf f}}^{{E_\varphi }}}$$ are calculated by formulae  \begin{equation*} G_{lj}^{{E_\varphi }} = \frac{1}{{{\mu _0}}}\left( {\int_{{{\Omega ^{rz}}}}{}\nabla {\psi _l} \cdot \nabla {\psi _j}r{\rm{d}}r{\rm{d}}z + \int_{{{\Omega ^{rz}}}}{}\frac{1}{r}{\psi _l}{\psi _j}{\rm{d}}r{\rm{d}}z} \right),\quad M_{lj}^{{E_\varphi }} = \int_{{{\Omega ^{rz}}}}{}\sigma _r^{p,rz}{\psi _l}{\psi _j}r{\rm{d}}r{\rm{d}}z, \end{equation*}   \begin{equation} f_l^{{E_\varphi }} = I{\psi _l}\left( {{r_\delta },{z_\delta }} \right){r_\delta }. \end{equation} (38)When FE approximation is used for this problem, we exclude a small subdomain in the vicinity of r = 0 from the computational domain Ωrz, on its boundary Γ0 the homogeneous Dirichlet boundary condition $${ {E_\varphi ^{}} |_{{\Gamma _0}}} = 0$$ is set (as on the distant boundaries of Ωrz). As mentioned above, for obtaining the field of the ‘vertical’ frame, it is necessary to solve two problems: eq. (20) and equation system (28)–(30). The equivalent variational formulation for eq. (20) takes the form   \begin{equation} \frac{1}{{{\mu _0}}}\int_{{{\Omega ^{rz}}}}{{\nabla A_\xi ^{} \cdot \nabla \nu r{\rm{d}}r{\rm{d}}z}} + i\omega \int_{{{\Omega ^{rz}}}}{{\sigma _r^{p,rz}A_\xi ^{}\nu r{\rm{d}}r{\rm{d}}z}} = Ih_L^{}\left( {\nu \left( {0,{z_1}} \right) - \nu \left( {0,{z_2}} \right)} \right), \end{equation} (39)where ν(0, z1) and ν(0, z2) are the results of delta-functions $${\delta ^{{z_1}}}$$ and $${\delta ^{{z_2}}}$$ operations on the test function ν(r, z). When this problem is solved in the vicinities of the points (0, z1) and (0, z2), FE solution can grow infinitely, and, therefore, when calculating $${{{\bf \vec{E}}}^{p,2}}$$, it is not allowed to come to these points closer than some pre-set distance. The solution $$A_\xi ^{}$$ of eq. (39) is also expressed as a linear combination of basis functions ψj(r, z)   \begin{equation} A_\xi ^{} = \sum\limits_{j = 1}^{{n^{rz}}} {q_j^{A_\xi ^{}}{\psi _j}} . \end{equation} (40)As a result of substituting (40) into (39) and replacing the test function with all basis functions ψl in turn, the following SLAE is obtained:   \begin{equation} \left( {{{{\bf G}}^{A_\xi ^{}}} + i\omega {{{\bf M}}^{A_\xi ^{}}}} \right){{{\bf q}}^{A_\xi ^{}}} = {{{\bf f}}^{A_\xi ^{}}}, \end{equation} (41)where the components of matrices $${{{\bf G}}^{A_\xi ^{}}}$$, $${{{\bf M}}^{A_\xi ^{}}}$$, and vector $${{{\bf f}}^{A_\xi ^{}}}$$ are calculated by formulae  \begin{equation*} G_{lj}^{A_\xi ^{}} = \frac{1}{{{\mu _0}}}\int_{{{\Omega ^{rz}}}}{}\nabla {\psi _l} \cdot \nabla {\psi _j}r{\rm{d}}r{\rm{d}}z,\quad M_{lj}^{A_\xi ^{}} = \int_{{{\Omega ^{rz}}}}{}\sigma _r^{p,rz}{\psi _l}{\psi _j}r{\rm{d}}r{\rm{d}}z, \end{equation*}   \begin{equation} f_l^{A_\xi ^{}} = Ih_L^{}\left( {{\psi _l}\left( {0,{z_1}} \right) - {\psi _l}\left( {0,{z_2}} \right)} \right). \end{equation} (42)The equivalent variational formulation for the equation system (28)–(30) takes the form   \begin{equation} \frac{1}{{{\mu _0}}}\left( {\int_{{{\Omega ^{rz}}}}{{\nabla \tilde{A}_r^{} \cdot \nabla \nu r{\rm{d}}r{\rm{d}}z}} + \int_{{{\Omega ^{rz}}}}{{\frac{1}{r}\tilde{A}_r^{}\nu {\rm{d}}r{\rm{d}}z}}} \right) + i\omega \int_{{{\Omega ^{rz}}}}{{\sigma _r^{p,rz}\tilde{A}_r^{}\nu r{\rm{d}}r{\rm{d}}z}} + \int_{{{\Omega ^{rz}}}}{{\sigma _r^{p,rz}\frac{{\partial V}}{{\partial r}}\nu r{\rm{d}}r{\rm{d}}z}} = 0, \end{equation} (43)  \begin{equation} \frac{1}{{{\mu _0}}}\int_{{{\Omega ^{rz}}}}{{\nabla A_z^{} \cdot \nabla \nu r{\rm{d}}r{\rm{d}}z}} + i\omega \int_{{{\Omega ^{rz}}}}{{\sigma _z^{p,rz}A_z^{}\nu r{\rm{d}}r{\rm{d}}z}} + \int_{{{\Omega ^{rz}}}}{{\sigma _z^{p,rz}\frac{{\partial V}}{{\partial z}}\nu r{\rm{d}}r{\rm{d}}z}} = \int_{{{\Omega ^{rz}}}}{{{J_z}\nu r{\rm{d}}r{\rm{d}}z}}, \end{equation} (44)  \begin{equation} \int_{{{\Omega ^{rz}}}}{{{\sigma ^{p,rz}}\nabla V \cdot \nabla \nu \,\,r{\rm{d}}r{\rm{d}}z}} + i\omega \int_{{{\Omega ^{rz}}}}{{{\sigma ^{p,rz}}{{( {\tilde{A}_r^{},A_z^{}} )}^{\rm{T}}} \cdot \nabla \nu \,\,r{\rm{d}}r{\rm{d}}z}} = i\omega \int_{{{\Omega ^{rz}}}}{{\sigma _r^{p,rz}u\nu r{\rm{d}}r{\rm{d}}z}}, \end{equation} (45)where ν(r, z) is the test function as before. Functions $$\tilde{A}_r^{}$$, $$A_z^{}$$ and V are represented in the form   \begin{equation} \tilde{A}_r^{} = \sum\limits_{i = 1}^{{n^{rz}}} {q_i^{\tilde{A}_r^{}}{\psi _i}} ,\quad A_z^{} = \sum\limits_{i = 1}^{{n^{rz}}} {q_i^{{A_z}}{\psi _i}} ,\quad V = \sum\limits_{i = 1}^{{n^{rz}}} {q_i^V{\psi _i}} . \end{equation} (46)Values $$q_i^{\tilde{A}_r^{}}$$, $$q_i^{{A_z}}$$ and $$q_i^V$$ are found from the SLAE which is obtained by substituting (46) into equation system (43)–(45) and replacing the test function ν with the basis function ψi in each equation. The obtained SLAE takes the following form:   \begin{equation} {{{\bf C}}^{AV}}{{{\bf q}}^{AV}} = {{{\bf f}}^{AV}}. \end{equation} (47)The size of system (47) is equal to 3nrz × 3nrz. The matrix CAV has a block structure with 3 × 3 block elements, vectors qAV and fAV also have a block structure with 3 × 1 block elements. The vector qAV has the following form: $${{{\bf q}}^{AV}} = {\{ {q_1^{\tilde{A}_r^{}},q_1^{A_z^{}},q_1^V,q_2^{\tilde{A}_r^{}},q_2^{A_z^{}},q_2^V,...} \}^{\rm{T}}}$$. The block elements of the matrix CAV take the following form:   \begin{equation} {{\bf C}}_{lj}^{AV} = \left[ {\begin{array}{c@{\quad}c@{\quad}c} {G_{ij}^{{{\tilde{A}}_r}} + i\omega M_{lj}^r}&0&{D_{lj}^r}\\ 0&{G_{ij}^{{A_z}} + i\omega M_{lj}^z}&{D_{lj}^z}\\ {i\omega D_{jl}^r}&{i\omega D_{jl}^z}&{G_{ij}^V} \end{array}} \right], \end{equation} (48)where   \begin{equation} G_{lj}^{{{\tilde{A}}_r}} = \frac{1}{{{\mu _0}}}\left( {\int_{{{\Omega ^{rz}}}}{}\nabla {\psi _l} \cdot \nabla {\psi _j}r{\rm{d}}r{\rm{d}}z + \int_{{{\Omega ^{rz}}}}{}\frac{1}{r}{\psi _l}{\psi _j}{\rm{d}}r{\rm{d}}z} \right),\quad G_{lj}^{{A_z}} = \frac{1}{{{\mu _0}}}\int_{{{\Omega ^{rz}}}}{}\nabla {\psi _l} \cdot \nabla {\psi _j}r{\rm{d}}r{\rm{d}}z, \end{equation} (49)  \begin{equation} G_{lj}^V = \int_{{{\Omega ^{rz}}}}{{{\sigma ^{p,rz}}\nabla {\psi _l} \cdot \nabla {\psi _j}\,\,r{\rm{d}}r{\rm{d}}z}},\quad D_{lj}^r = \int_{{{\Omega ^{rz}}}}{}\sigma _r^{p,rz}{\psi _l}\frac{{\partial {\psi _j}}}{{\partial r}}r{\rm{d}}r{\rm{d}}z,\quad D_{lj}^z = \int_{{{\Omega ^{rz}}}}{}\sigma _z^{p,rz}{\psi _l}\frac{{\partial {\psi _j}}}{{\partial z}}r{\rm{d}}r{\rm{d}}z, \end{equation} (50)  \begin{equation} M_{lj}^r = \int_{{{\Omega ^{rz}}}}{}\sigma _r^{p,rz}{\psi _l}{\psi _j}r{\rm{d}}r{\rm{d}}z,\quad M_{lj}^z = \int_{{{\Omega ^{rz}}}}{}\sigma _z^{p,rz}{\psi _l}{\psi _j}r{\rm{d}}r{\rm{d}}z. \end{equation} (51)The block elements of the vector fAV take the following form:   \begin{equation} {{\bf f}}_l^{AV} = \left[ {\begin{array}{@{}*{1}{c}@{}} 0\\ {f_l^z}\\ {i\omega f_l^u} \end{array}} \right],\quad f_l^z = \int_{{{\Omega ^{rz}}}}{{{J_z}{\psi _l}r{\rm{d}}r{\rm{d}}z}},\quad f_l^u = \int_{{{\Omega ^{rz}}}}{{\sigma _r^{p,rz}u{\psi _l}r{\rm{d}}r{\rm{d}}z}}. \end{equation} (52)Note that when calculating $$f_l^u$$ in (52), the function u is found as $$u = \sum_{j = 1}^{{n^{rz}}} {q_j^u{\psi _j}} $$ by solving the equation, which is analogous to eq. (20) (or 39), that is, with the use of eqs (41) and (42) when $$h_L^{} = 1$$. Therefore, the component $$f_l^u$$ in (52) can be calculated as $$f_l^u = \sum_j {M_{lj}^rq_j^u} $$ where $$M_{lj}^r$$ is defined by (51). There is a further point to be made here. In order to obtain the complex symmetrical SLAE, every third row of the system must be divided by iω. 3.2 Finite-element approximation for mathematical models used for calculating the 3-D field in anisotropic medium The equivalent variational formulation for eq. (1) takes the following form:   \begin{eqnarray} \int_{\Omega }{{\frac{1}{{{\mu _0}}}}}( {\nabla \times {{{{\bf \vec{E}}}}^s}} ) \cdot ( {\nabla \times {{\bf \vec{\boldsymbol\psi }}}} ){\rm{d}}\Omega + i\omega \int_{\Omega }{{\sigma {{{{\bf \vec{E}}}}^s}}} \cdot {{\bf \vec{\boldsymbol\psi }}}{\rm{d}}\Omega - {\omega ^2}\int_{\Omega }{{\varepsilon {{{{\bf \vec{E}}}}^s}}} \cdot {{\bf \vec{\boldsymbol\psi }}}{\rm{d}}\Omega &=& - i\omega \int_{\Omega }{{\left( {\sigma - {\sigma ^p}} \right){{{{\bf \vec{E}}}}^p}}} \cdot {{\bf \vec{\boldsymbol\psi }}}{\rm{d}}\Omega + {\omega ^2}\int_{\Omega }{{\varepsilon {{{{\bf \vec{E}}}}^p}}} \cdot {{\bf \vec{\boldsymbol\psi }}}{\rm{d}}\Omega \,\,, \end{eqnarray} (53)where Ω is the computational domain (and dΩ = dxdydz), the test vector function $${{\bf \vec{\boldsymbol\psi }}}$$ and the unknown vector function $${{{\bf \vec{E}}}^s}$$ are the members of the space Hcurl(Ω) of the vector functions $${{\bf \vec{\Phi }}}$$ with continuous tangential components, which satisfy conditions $$\int_{\Omega }{{{{\bf \vec{\Phi }}} \cdot {{\bf \vec{\Phi }}}{\rm{d}}\Omega }} < \infty $$ and $$\int_{\Omega }{{( {\nabla \times {{\bf \vec{\Phi }}}} ) \cdot ( {\nabla \times {{\bf \vec{\Phi }}}} ){\rm{d}}\Omega }} < \infty $$. The function $${{{\bf \vec{E}}}^s}$$ is expressed as   \begin{equation} {{{\bf \vec{E}}}^s} = \sum\limits_{j = 1}^n {q_j^s{{\vec{\boldsymbol\psi }}_j}} , \end{equation} (54)where $${\vec{\boldsymbol\psi }_j}$$ are the basis vector functions, n is their number (which coincides with the number of the edges in the 3-D mesh). As a result of substituting (54) into (53) and replacing the test function $${{\bf \vec{\boldsymbol\psi }}}$$ with the basis functions $${\vec{\boldsymbol\psi }_l}$$ in turn, we obtain the following SLAE:   \begin{equation} \left( {{{\bf G}} + i\omega {{{\bf M}}^\sigma } - {\omega ^2}{{{\bf M}}^\varepsilon }} \right){{{\bf q}}^s} = {{\bf f}}, \end{equation} (55)where the components of matrices G, Mσ and Mε, and the vector f are defined by relations   \begin{equation} {G_{lj}} = \frac{1}{{{\mu _0}}}\int_{\Omega }{}\left( {\nabla \times {{\vec{\boldsymbol\psi }}_l}} \right) \cdot \left( {\nabla \times {{\vec{\boldsymbol\psi }}_j}} \right){\rm{d}}\Omega ,\quad M_{lj}^\sigma = \int_{\Omega }{}\sigma {\vec{\boldsymbol\psi }_l} \cdot {\vec{\boldsymbol\psi }_j}{\rm{d}}\Omega ,\quad M_{lj}^\varepsilon = \int_{\Omega }{}\varepsilon {\vec{\boldsymbol\psi }_l} \cdot {\vec{\boldsymbol\psi }_j}{\rm{d}}\Omega , \end{equation} (56)  \begin{equation} {f_l} = - i\omega \int_{\Omega }{}\left( {\sigma - {\sigma ^p}} \right){{{\bf \vec{E}}}^p} \cdot {\vec{\boldsymbol\psi }_l}{\rm{d}}\Omega + {\omega ^2}\int_{\Omega }{{\varepsilon {{{{\bf \vec{E}}}}^p}}} \cdot {\vec{\boldsymbol\psi }_l}{\rm{d}}\Omega . \end{equation} (57)The function $${{{\bf \vec{E}}}^p}$$ is expressed as a linear combination of the basis functions $${\vec{\boldsymbol\psi }_j}$$  \begin{equation} {{{\bf \vec{E}}}^p} = \sum\limits_{j = 1}^n {q_j^p{{\vec{\boldsymbol\psi }}_j}} , \end{equation} (58)where expansion coefficients $$q_j^p$$ in basis $$\{ {{{\vec{\boldsymbol\psi }}_j}} \}$$ are the values of tangential components of $${{{\bf \vec{E}}}^p}$$ on the FEs’ edges. In this case, the vector f in the right-hand side of FE equation system (55) can be calculated by the formula   \begin{equation} {{\bf f}} = \left( { - i\omega {{{{\bf \bar{M}}}}^\sigma } + {\omega ^2}{{{\bf M}}^\varepsilon }} \right){{\bf q}}_{}^p, \end{equation} (59)where components of the matrix $${{{\bf \bar{M}}}^\sigma }$$ are defined by the formula   \begin{equation} \bar{M}_{lj}^\sigma = \int_{\Omega }{}\left( {\sigma - {\sigma ^p}} \right){\vec{\boldsymbol\psi }_l} \cdot {\vec{\boldsymbol\psi }_j}{\rm{d}}\Omega . \end{equation} (60) 4 SPECIAL ASPECTS OF IMPLEMENTATION For the solution of axisymmetrical problems (when we calculate the primary field $${{{\bf \vec{E}}}^p}$$), we use the regular rectangular meshes with smoothly increasing of the step when moving away from the source, which is located on the z-axis. The step in the vicinity of the source is set about 0.01 m, the factor χ of increase of the mesh step when moving away from the source in the coordinates R and Z is set in the range of 1.05–1.08 (i.e. the value of the next mesh step hnext in each of the coordinates is defined by the value of the previous step hprev in the corresponding coordinate as hnext = χ ⋅ hprev). The approach proposed in the previous sections for modelling the primary field requires much lower computational costs in comparison with computational costs required for calculating the 3-D secondary field. Therefore, we use the special algorithms of optimization only for the 3-D meshes. These algorithms make it possible to reduce the computational time required for calculating the secondary field by approximately an order in comparison with the calculations on regular meshes (which are generated according to the principle considered above for 2-D meshes). Therefore, in this Section, we consider the algorithms of generating the optimized irregular non-conforming 3-D meshes and principles of constructing the FE approximation on them, as well as another aspects (grouping of logging tool positions and automatic choice of the background medium), which allow additional significant reduction of the computational costs. Note that we do not use the approaches with adaptive mesh generation, which allow the convergence towards the solution using the sequence of meshes. These approaches are based on various error estimators (Schwarzbach et al.2011; Grayver & Kolev 2015; Yin et al.2016, etc.) and require the additional costs caused by the calculations on several meshes in order to obtain the adaptive mesh and for each transmitter position. At the same time, these approaches do not always guarantee the required accuracy when the signals are calculated in the receivers. This fact is implicitly confirmed by Grayver & Kolev (2015) where, after generating the adaptive mesh (and satisfying the corresponding criteria), the additional refinements are performed in the vicinities of the receivers. 4.1 Generation of the optimized finite-element meshes For the FE solution, we use the optimized meshes with local refinements. Note that the trajectory length in the logging problems is usually much bigger than the size of the computational domain, in which the field for one position of the source can be computed with the required accuracy. For increasing the numerical solution accuracy it is necessary to refine the mesh in the vicinity of the transmitters and receivers. Therefore, during the calculation of the signals received by the logging tool when it is moving along the trajectory, it is necessary to change not only the positions of the local refinements of the mesh but also the computational domain itself (the analogous approach has been considered in (Yang et al.2014), but it uses regular meshes). The accuracy of approximation of the computational domain geometry is the important aspect for obtaining the high accuracy solution. On the one hand, this approximation should be quite good in the vicinity of the logging tool (i.e. the error caused by the inaccurate boundaries approximation should not exceed the error of the numerical solution approximation). On the other hand, the excessive mesh refinement caused by good approximation of the boundaries, which are distant from the logging tool, can lead to the significant increase in the number of FEs and, consequently, computational costs. The numerical results have shown that high accuracy geometry approximation is required only in some subdomain around the logging tool, and rough geometry approximation can be used at some distance from it. Note that though we use the irregular non-conforming meshes, we do not use FEs with curvilinear or inclined boundaries. This is related to the fact that, if we use the FEs with curvilinear or inclined boundaries, which smoothly approximate the corresponding boundaries of the domain, the application of the primary-secondary field approach can lead to the intersection of FEs of the 3-D mesh by the layers interfaces of the background medium. In this case, the secondary field $${{{\bf \vec{E}}}^s}$$ should be discontinuous inside those FEs which are intersected by the interfaces of background medium layers with different conductivity (because $${{{\bf \vec{E}}}^t} = {{{\bf \vec{E}}}^p} + {{{\bf \vec{E}}}^s}$$ inside the FE should be continuous, but $${{{\bf \vec{E}}}^p}$$ will be discontinuous if the FE is intersected by the layers with a different conductivity). Therefore, we use another approach (stair-like approximation) in which the FEs are never intersected by the layers interfaces of the background medium. Mesh generation algorithm consists of four main stages. At the first stage, the size and structure of the computational domain are defined depending on the logging tool position (or several positions; it will be described in Section 4.3). The second stage is the stage of generating the so-called base mesh, which does not include inner borders of the computational domain and is refined in the subdomain containing the logging tool. At the third stage, the geoelectrical inhomogeneities are embedded into the mesh taking into account the requirements to the accuracy of their boundaries approximation depending on the distance from the logging tool. And, finally, at the fourth stage, the local refinements in the vicinity of receivers and transmitters are performed. Thus, when generating the base mesh, we select the parallelepiped subdomain, which contains the transmitters and receivers locations for one tool position (or several tool positions included in one group). Inside this subdomain, the base mesh is generated with a constant step in all directions, and outside this subdomain, the cells size increases smoothly with moving away from the selected subdomain boundaries. Thus, the base mesh for a group of tool positions on the trajectory is defined by three parameters. They are the following: the step inside the selected subdomain, the distance from the tool to boundaries of the computational domain, and the factor of increase of the cell size (when moving away from the selected subdomain boundaries). Based on multiple numerical experiments performed for different geoelectrical models, we have obtained the following values of these parameters. For the frequency of 100 kHz, which we take as a ‘basic’ one, the size of the minimal step in the base mesh is set equal to 0.5 m, and the distance from the tool to the boundaries of the computational domain is set equal to 20 m in all directions. For other frequencies, these sizes are scaled by multiplying by the square root of the frequencies ratio. The factor of increase of the cells size (when moving away from the boundaries of the subdomain containing all locations of transmitters and receivers included in a group) is taken equal to 1.2. Before the third stage of the algorithm, all geoelectrical inhomogeneities included in the computational domain (or the whole computational domain) are approximated by parallelepipeds the size of which is defined by the step of geometry approximation. The step defining the approximation of geometry of embedded geoelectrical inhomogeneities in the vicinity of receivers’ locations is about 0.2 m for the problems without a borehole and invasion zone, and it is about 0.05 m for the problems in which the borehole and invasion zone need to be taken into account. Outside the vicinity of receivers’ locations, the step of boundaries approximation increases proportionally to the increase in the step size of the base mesh. The local refinements in the vicinity of receivers (at the fourth stage of the algorithm) are defined by the following parameters: the number of refinements inside the FEs into which the receivers get, the number of neighbouring levels of FEs, and the number of refinements for each neighbouring levels. For the examples considered in this paper, the values of these parameters are taken equal to 2, 3 and 1, respectively. As mentioned above, for solving the problems considered here, we use the mathematical models based on the representation of the total field as a sum of the primary (corresponding to the field induced by the source in the horizontally layered background medium) and secondary (which is the impact field of 3-D and 2-D geoelectrical inhomogeneities) components. Note that the local refinements in the vicinity of the transmitter should be performed only in the case if the transmitter is located in the subdomain (of the computational domain) in which σp ≠ σ. As an example, consider the geoelectrical model the vertical section of which is shown in Fig. 3. The logging tool trajectory is shown by black solid line. Figure 3. View largeDownload slide The geoelectrical model specific for induction logging problems. The model contains five layers (interfaces of which are slightly sloped) and the fault. The layers of equivalent conductivity are shown by the same colour. The conductivity of the middle layer is anisotropic: lateral conductivity is 0.256 Sm m−1 and vertical conductivity is 0.025 Sm m−1. Figure 3. View largeDownload slide The geoelectrical model specific for induction logging problems. The model contains five layers (interfaces of which are slightly sloped) and the fault. The layers of equivalent conductivity are shown by the same colour. The conductivity of the middle layer is anisotropic: lateral conductivity is 0.256 Sm m−1 and vertical conductivity is 0.025 Sm m−1. We consider the logging tool that includes one transmitter and two receivers the offsets of which are 4.5 and 9.5 m. The distance between two adjacent logging tool positions is 1 m. The value of the frequency is equal to 100 kHz. Fig. 4(a) shows the base mesh, which is obtained after the first and second stages of the algorithm, for one group of three neighbouring logging tool positions on the trajectory. Fig. 4(b) shows the mesh obtained after the third stage of the algorithm, that is, after embedding 2-D geoelectrical inhomogeneities into the mesh. Fig. 5(a) shows the mesh obtained after the fourth stage of the algorithm, that is, after local refinement in the vicinity of the receivers and transmitters. Fig. 5(b) shows the 3-D view of the sectioned FE mesh. In Fig. 6, the FE mesh is shown for other logging tool positions on the trajectory. Figure 4. View largeDownload slide (a) The vertical section of the base mesh obtained after the first and second algorithm stages and (b) the vertical section of the mesh obtained after embedding 2-D geoelectrical inhomogeneities into the mesh. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the ith transmitter. Figure 4. View largeDownload slide (a) The vertical section of the base mesh obtained after the first and second algorithm stages and (b) the vertical section of the mesh obtained after embedding 2-D geoelectrical inhomogeneities into the mesh. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the ith transmitter. Figure 5. View largeDownload slide (a) The vertical section of the mesh obtained after local refinements in the vicinities of receivers and transmitters and (b) 3-D view of the sectioned finite-element mesh. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the ith transmitter. Figure 5. View largeDownload slide (a) The vertical section of the mesh obtained after local refinements in the vicinities of receivers and transmitters and (b) 3-D view of the sectioned finite-element mesh. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the ith transmitter. Figure 6. View largeDownload slide The vertical sections of the meshes generated for different groups of the adjacent tool positions. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the jth transmitter. Figure 6. View largeDownload slide The vertical sections of the meshes generated for different groups of the adjacent tool positions. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the jth transmitter. 4.2 Finite-element approximation on the meshes with hanging nodes 4.2.1 ‘Non-conforming’ and ‘conforming’ bases on the meshes with hanging nodes The meshes presented in the previous section contain the non-conforming FEs. On these meshes, we introduce two sets of the basis functions: $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$ and $$\{ {\vec{\boldsymbol\psi }_k^c} \}$$, where the superscript ‘nc’ denotes non-conforming, and superscript ‘c’ denotes conforming. The basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$ are associated with all edges of the mesh with hanging nodes. On each FE Ωe containing the corresponding (kth) edge, $$\vec{\boldsymbol\psi }_k^{nc}$$ is defined in the standard way: at this edge, the basis function $$\vec{\boldsymbol\psi }_k^{nc}$$ has the only non-zero (equal to 1) component, which corresponds to the coordinate (x, y or z) coinciding with the edge direction, and, on Ωe, this component decreases to zero as a bilinear function of two other coordinates when coming to three other edges with the same direction. Obviously, the functions $$\vec{\boldsymbol\psi }_k^{nc} \notin {{{\bf H}}^{{\rm{curl}}}}( \Omega )$$ occur in the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. As an example, consider the fragment of the irregular non-conforming mesh (with hanging nodes), which is shown in Fig. 7. Denote the basis function, which is associated with the edge connecting nodes with numbers 5 and 26, as $$\vec{\boldsymbol\psi }_5^{nc}$$. This function is non-zero on the FEs Ω1, Ω2, and Ω3, and zero on all other FEs. Further, we will denote the edges, FEs, and their faces by the sets of nodes, which are listed in braces. Thus, the basis function $$\vec{\boldsymbol\psi }_5^{nc}$$ is associated with the edge {5, 26}, FE Ω1 = {1, 2, 4, 5, 22, 23, 25, 26}, FE Ω1 has the faces {1, 2, 4, 5}, {1, 2, 22, 23}, etc. As an example, consider the basis functions $$\vec{\boldsymbol\psi }_{36}^{nc}$$, $$\vec{\boldsymbol\psi }_{37}^{nc}$$ and $$\vec{\boldsymbol\psi }_{38}^{nc}$$, which are associated with the edges {5, 10}, {10, 16} and {16, 26}, respectively. These edges are the sub-edges of the edge {5, 26}. The basis function $$\vec{\boldsymbol\psi }_{36}^{nc}$$ is non-zero only on FE Ω4 = {5, 6, 8, 9, 10, 12, 13, 15}, basis function $$\vec{\boldsymbol\psi }_{37}^{nc}$$ is non-zero only on FE Ω5 = {10, 11, 13, 14, 16, 17, 19, 20}, basis function $$\vec{\boldsymbol\psi }_{38}^{nc}$$ is non-zero only Ω7 = {16, 17, 19, 20, 26, 27, 30, 31}, and on all other FEs, these basis functions are equal to zero. Figure 7. View largeDownload slide The fragment of the irregular non-conforming finite-element mesh with hanging nodes. Figure 7. View largeDownload slide The fragment of the irregular non-conforming finite-element mesh with hanging nodes. Besides the basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$ associated with the sub-edges of ‘big’ edges (which are the edges of other FEs), in the ‘non-conforming’ basis, there are functions, which are associated with the FEs’ edges lying within the faces of other FEs (naturally, such edges do not coincide with the boundaries of these faces). In Fig. 7, these are the functions $$\vec{\boldsymbol\psi }_{39}^{nc}$$ and $$\vec{\boldsymbol\psi }_{40}^{nc}$$, which are associated with the edges {11, 17} and {17, 27}, respectively. The function $$\vec{\boldsymbol\psi }_{39}^{nc}$$ is non-zero only on FE Ω5 (see above its nodes) and Ω6 = {11, 12, 14, 15, 17, 18, 20, 21}, and the function $$\vec{\boldsymbol\psi }_{40}^{nc}$$ is non-zero only on FE Ω7 and Ω8 = {17, 18, 20, 21, 27, 28, 31, 32}. In general, by the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$ we mean the basis functions associated with all edges (including sub-edges). These basis functions are defined in the standard form on the FEs including these edges. On the non-conforming FE mesh with hanging nodes, we introduce the basis functions $$\vec{\boldsymbol\psi }_j^c$$ with continuous tangential components in the computational domain Ω, that is, $$\vec{\boldsymbol\psi }_j^c \in {{{\bf H}}^{{\rm{curl}}}}( \Omega )$$. These basis functions are associated only with the ‘main’ edges of the mesh. By the ‘main’ edges we mean such edges that are not the sub-edges and do not lie within the faces of any FEs. Define the basis functions $$\vec{\boldsymbol\psi }_j^c$$ in the form of linear combination of the basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$  \begin{equation} \vec{\boldsymbol\psi }_j^c = \sum\limits_k {{T_{jk}}\vec{\boldsymbol\psi }_k^{nc}} , \end{equation} (61)that is, the basis $$\{ {\vec{\boldsymbol\psi }_j^c} \}$$ is constructed from the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$ with the help of the transition matrix T = {Tjk}. The size of this matrix is nc × n, where n is the total number of the edges in the mesh (including sub-edges and edges, which lie within the faces of any FEs), and nc is the number of ‘main’ edges only. Thus, specifically with the help of the transition matrix T, we construct the basis $$\{ {\vec{\boldsymbol\psi }_j^c} \}$$ in which all vector functions have continuous tangential components in Ω (i.e. ∀j$$\vec{\boldsymbol\psi }_j^c \in {{{\bf H}}^{{\rm{curl}}}}( \Omega )$$). We note a very important aspect related to the ‘non-conforming’ basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. Besides the basis functions with the discontinuous tangential components (which we have called ‘non-conforming’), the ‘conforming’ basis functions (with continuous tangential components) are also included into the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. They are associated with the edges, which we will name ‘regular’. These edges do not contain subedges, are not subedges themselves, and do not lie within any face. For example, in Fig. 7, the ‘regular’ edges are the edges {1, 22}, {4, 25}, {2, 23}, etc. (if we assume that either the faces of FE Ω1 lie on the computational domain boundary, or FE Ω1 and FEs adjacent to Ω1, which are not shown in this mesh fragment, are conforming). The superscript ‘nc’ added to these basis functions (associated with ‘regular’ edges) denotes only that they are included in the basis in which there are vector functions with discontinuous tangential components. Thus, the superscript ‘nc’ added to the function from the ‘non-conforming’ basis does not always denote that this function necessarily have discontinuous tangential components, that is, among the functions from $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$, there are functions with continuous tangential components, which are associated with the ‘regular’ edges. Therefore, further, the superscript ‘nc’ added to the basis function $$\vec{\boldsymbol\psi }_k^{nc}$$ will not rigorously denote that this function necessarily has discontinuous tangential components (i.e. the term ‘conforming’ function $$\vec{\boldsymbol\psi }_k^{nc}$$ is applicable and denotes that such a function from the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$ is associated with the ‘regular’ edge). However, all ‘regular’ edges are the ‘main’, but unlike ‘regular’, the ‘main’ edges can contain the sub-edges (which are the edges of some FEs). For example, in Fig. 7, the edge {5, 26} is ‘main’, but not ‘regular’, because it contains subedges {5, 10} (the edge of FE Ω4), {10, 16} (the edge of FE Ω5), and {16, 26} (the edge of FE Ω7). 4.2.2 The construction of the ‘conforming’ basis from the ‘non-conforming’ Consider how the ‘conforming’ basis function $$\vec{\boldsymbol\psi }_j^c$$ can be constructed from the ‘non-conforming’ functions $$\vec{\boldsymbol\psi }_k^{nc}$$. In so doing, we really develop the algorithm of generating the matrix T. Let the face of one FE (this face will be called ‘big’) be adjacent to the several faces of other FEs (these faces will be called ‘small’). Then, in order to ensure on this ‘big’ face the continuity of the tangential component of the basis function $$\vec{\boldsymbol\psi }_j^c$$ associated with one of its edges, it is necessary to define the weights Tjk in the expansion (61) in basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$. The weights Tjk are equal to the values of the tangential component $$\vec{\boldsymbol\psi }_j^c$$ on the corresponding edges of the ‘small’ faces. We illustrate this by the example of the mesh shown in Fig. 7. Consider the face {5, 6, 26, 28} of FE Ω2 (this is the ‘big’ face) adjacent to the faces {5, 6, 10, 12} of FE Ω4, {10, 11, 16, 17} of FE Ω5, {11, 12, 17, 18} of FE Ω6, {16, 17, 26, 27} of FE Ω7 and {17, 18, 27, 28} of FE Ω8 (these are the ‘small’ faces). Let, as before, the vector function associated with the edge {5, 26} have number 5 in the ‘non-conforming’ basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$, and vector functions associated with the edges {5, 10}, {10, 16}, {16, 26}, {11, 17} and {17, 27} have numbers 36, 37, 38, 39 and 40, respectively. In the ‘conforming’ basis $$\{ {\vec{\boldsymbol\psi }_j^c} \}$$, the vector function associated with the edge {5, 26} has number 5 too. Then, the ‘conforming’ vector function associated with the edge {5, 26} takes the form   \begin{equation} \vec{\boldsymbol\psi }_5^c = \vec{\boldsymbol\psi }_5^{nc} + \vec{\boldsymbol\psi }_{36}^{nc} + \vec{\boldsymbol\psi }_{37}^{nc} + \vec{\boldsymbol\psi }_{38}^{nc} + \frac{{{h_{11,12}}}}{{{h_{10,12}}}}\left( {\vec{\boldsymbol\psi }_{39}^{nc} + \vec{\boldsymbol\psi }_{40}^{nc}} \right), \end{equation} (62)where hν, η is the length of the edge {ν,η}. Coefficients T5, 5 = 1, T5, 36 = 1, T5, 37 = 1, T5, 38 = 1, $${T_{5,39}} = {{{h_{11,12}}} / {{h_{10,12}}}}$$, and $${T_{5,40}} = {{{h_{11,12}}} / {{h_{10,12}}}}$$ are the values of the function $$\vec{\boldsymbol\psi }_5^c$$ at the edges {5, 10}, {10, 16}, {16, 26}, {11, 17} and {17, 27}, which lie on the ‘big’ face {5, 6, 26, 28} and with which the functions $$\vec{\boldsymbol\psi }_{36}^{nc}$$, $$\vec{\boldsymbol\psi }_{37}^{nc}$$, $$\vec{\boldsymbol\psi }_{38}^{nc}$$, $$\vec{\boldsymbol\psi }_{39}^{nc}$$ and $$\vec{\boldsymbol\psi }_{40}^{nc}$$, are associated. Note that the function $$\vec{\boldsymbol\psi }_5^c$$ identically coincides with the function $$\vec{\boldsymbol\psi }_5^{nc}$$ on the ‘big’ face, because on FE Ω2, both $$\vec{\boldsymbol\psi }_5^c$$ and $$\vec{\boldsymbol\psi }_5^{nc}$$ are defined by the same shape function. We also note that formula (62) defines the function $$\vec{\boldsymbol\psi }_5^c$$ only inside the FEs, and its values on the faces are defined as continuous extensions of its values from the FEs. On the face {5, 6, 26, 28} itself, formula (62) should be interpreted only in the sense of continuous extensions of the basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$ onto this face from the FEs, which contain the edges corresponding to $$\vec{\boldsymbol\psi }_k^{nc}$$, because all functions $$\vec{\boldsymbol\psi }_k^{nc}$$ in the right-hand side of (62) are discontinuous on this face. As a result, we really obtain the function $$\vec{\boldsymbol\psi }_5^c$$ with continuous tangential components on the face {5, 6, 26, 28}. Thus, we propose the following algorithm of calculating the matrix T, which makes it possible to obtain the ‘conforming’ basis functions $$\vec{\boldsymbol\psi }_j^c$$ from the ‘non-conforming’ (but standard for the vector FEM) basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$. To begin with, we assume that in the ‘non-conforming’ basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$, the functions associated with the ‘main’ edges are numbered first. The same numbering of the ‘main’ edges is used for the functions from the ‘conforming’ basis $$\{ {\vec{\boldsymbol\psi }_j^c} \}$$. Then, the left block of the matrix T with the size nc × nc is the identity matrix (i.e. Tjj = 1, and Tjk = 0 if k ≠ j). This is convenient both for the storage of the matrix T (i.e. in fact, there is no need to store this block) and for assembling the matrices of FEs. Note that the basis functions $$\vec{\boldsymbol\psi }_j^c$$ associated with the ‘regular’ edges are identically equal to $$\vec{\boldsymbol\psi }_j^{nc}$$, and for these edges, all components of the corresponding (jth) rows in the right block of the matrix T are equal to zero (i.e. Tjk = 0 if k > nc). Non-zero components Tjk if k > nc are calculated in the following way. If the edge numbered k is the sub-edge of the edge numbered α ≤ nc, then Tαk = 1, and all other components of the kth column of the matrix T are zero (i.e. ∀β ≠ α Tβk = 0). If the edge numbered k lies within the ‘big’ face, on which the basis functions $$\vec{\boldsymbol\psi }_\alpha ^{nc}$$ and $$\vec{\boldsymbol\psi }_\beta ^{nc}$$ for α, β ≤ nc have non-zero component tangential to this edge, then only two components Tαk and Tβk are non-zero in the kth column of the matrix T, and their values are equal to the values of functions $$\vec{\boldsymbol\psi }_\alpha ^{nc}$$ and $$\vec{\boldsymbol\psi }_\beta ^{nc}$$ at the kth edge, respectively. Thus, for the mesh shown in Fig. 7, if the edge {11, 17}, which lies within the face {5, 6, 26, 28}, has number 39, and edges {5, 26} and {6, 28} have numbers 5 and 6, respectively, then only the components $${T_{5,39}} = {{{h_{11,12}}} / {{h_{10,12}}}}$$ and $${T_{6,39}} = {{{h_{10,11}}} / {{h_{10,12}}}}$$ (where, as before, hν, η is the length of the edge {ν,η}) are non-zero in the 39th column of the matrix T. Finally, we consider the situation when the edge numbered k > nc lies within the ‘big’ face and even if one of the edges of this ‘big’ face, which are parallel to the k-th edge, is not the ‘main’, that is, its number is α > nc. This situation is shown in Fig. 8. The edge {10, 16} lies within the face {4, 5, 24, 26} and, therefore, this face is ‘big’ over the edge {10, 16} and FEs {9, 10, 12, 13, 15, 16, 18, 19} and {10, 11, 13, 14, 16, 17, 19, 20}, which contain it. But the edge {4, 24} of this ‘big’ face lies within the face {2, 7, 22, 28}, which is ‘big’ over this edge and FE {2, 3, 4, 5, 22, 23, 24, 26}, which contains this edge. Therefore, the edge {4, 24} is not ‘main’ and its number is α > nc. In this case, the non-zero components of the kth column of the matrix T are calculated with the use of a ‘chain’. First, the values $${\tilde{T}_{\alpha k}}$$ are calculated as values of tangential components of the basis function $$\vec{\boldsymbol\psi }_\alpha ^{nc}$$ at the kth edge. For the example shown in Fig. 8 and edges {10, 16} (its number is equal to k) and {4, 24} (its number is equal to α), $${\tilde{T}_{\alpha k}}$$ is equal to $${{{h_{10,11}}} / {{h_{9,11}}}}$$. The tilde ‘∼’ over Tαk have been added because if α > nc, there are no such components in the matrix T (because its size is nc × n), and the values $${\tilde{T}_{\alpha k}}$$ if α > nc are really fictitious and used only in the ‘chain’ for calculating the non-fictitious components of the matrix T (i.e. components Tjk if j ≤ nc). Figure 8. View largeDownload slide The fragment of the irregular non-conforming finite-element mesh with hanging nodes such that the ‘chains’ appear when the transition matrix is generated. Figure 8. View largeDownload slide The fragment of the irregular non-conforming finite-element mesh with hanging nodes such that the ‘chains’ appear when the transition matrix is generated. Further, the ‘big’ face within which the edge numbered α lies is processed. For the edges of this face, which are parallel to the edge numbered α, the values of basis functions $$\vec{\boldsymbol\psi }_\gamma ^{nc}$$ and $$\vec{\boldsymbol\psi }_\nu ^{nc}$$ associated with these edges are calculated at the edge numbered α. If γ ≤ nc and ν ≤ nc, that is, the edges numbered γ and ν are the ‘main’ ones, then the components $${T_{\gamma k}} = {T_{\gamma \alpha }} \cdot {\tilde{T}_{\alpha k}}$$ and $${T_{\nu k}} = {T_{\nu \alpha }} \cdot {\tilde{T}_{\alpha k}}$$ are non-zero components in the kth column of the matrix T. Thus, for the mesh shown in Fig. 8, if the edges {10, 16}, {2, 22} and {7, 28} are numbered k, γ and ν, respectively, then $${T_{\gamma ,k}} = \frac{{{h_{10,11}}}}{{{h_{9,11}}}} \cdot \frac{{{h_{4,7}}}}{{{h_{2,7}}}}$$ and $${T_{\nu ,k}} = \frac{{{h_{10,11}}}}{{{h_{9,11}}}} \cdot \frac{{{h_{2,4}}}}{{{h_{2,7}}}}$$. In so doing, the generation of the kth column of the matrix T is completed. If at least one of the edges is not the ‘main’ one (i.e. γ > nc and/or ν > nc), then for this edge (or these two edges), instead of Tγα and/or Tνα the ‘fictitious’ values $${\tilde{T}_{\gamma \alpha }}$$ and/or $${\tilde{T}_{\nu \alpha }}$$ are calculated, and the corresponding ‘chains’ (or one of them) must be continued in a similar manner. Thus, the generation of the kth column of the matrix T will be completed when all ‘chains’, which have already been begun, are finished (i.e. the ‘main’ edge should be at the end of each ‘chain’). Note that a number of non-zero components of the matrix T (even without taking into account the identity submatrix with the size of nc × nc) is small. Therefore, for its storage it is reasonable to use the sparse format as for the global matrix of the FE equations system. The use of the matrix T allows sufficiently convenient generation of matrices G, Mσ, and Mε of FE equations system (55) for the ‘conforming’ basis $$\{ {\vec{\boldsymbol\psi }_j^c} \}$$ by assembling the standard FE matrices corresponding to the ‘non-conforming’ basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. 4.2.3 Generation of global matrices of finite-element equations system We find the solution $${{{\bf \vec{E}}}^s}$$ of eq. (53) as a linear combination of basis functions $$\vec{\boldsymbol\psi }_j^c$$:  \begin{equation*} {{{\bf \vec{E}}}^s} = \sum\limits_{j = 1}^{{n_c}} {q_j^s\vec{\boldsymbol\psi }_j^c} . \end{equation*} Note once again that basis functions $$\vec{\boldsymbol\psi }_j^c$$ have continuous tangential components in the whole computational domain Ω ($$\vec{\boldsymbol\psi }_j^c \in {{{\bf H}}^{{\rm{curl}}}}( \Omega )$$) and, therefore, are suitable for variational equation (53). We show how the matrices G, Mσ, and Mε, which define the matrix of FE equations system (55), can be calculated using the FE matrices, the components of which are obtained with the use of the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. The basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$ is defined by the standard shape vector functions on FEs, including the ‘non-conforming’ ones, that is, on each FE the basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$, which are non-zero on it, are the shape functions. Thus, shape functions are the local basis functions (on FE) of the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. We obtain the expressions for calculating the components of the matrix G (see eq. 56) with the use of the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. Considering that tangential components of basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$ can be discontinuous on the faces of the FEs Ωe, before turning to the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$, we represent the integral over $$\Omega = \bigcup\limits_e {{\Omega _e}} $$ in (56) as a sum of integrals over FEs Ωe, and after that we use expression (61):  \begin{eqnarray*} {G_{lj}} &=& \frac{1}{{{\mu _0}}}\int_{\Omega }{}\left( {\nabla \times \vec{\boldsymbol\psi }_l^c} \right) \cdot \left( {\nabla \times \vec{\boldsymbol\psi }_j^c} \right){\rm{d}}\Omega = \frac{1}{{{\mu _0}}}\sum\limits_e {\int_{{{\Omega _e}}}{}\left( {\nabla \times \vec{\boldsymbol\psi }_l^c} \right) \cdot \left( {\nabla \times \vec{\boldsymbol\psi }_j^c} \right){\rm{d}}\Omega }\nonumber \\ &=& \frac{1}{{{\mu _0}}}\sum\limits_e {\int_{{{\Omega _e}}}{}\left( {\nabla \times \sum\limits_k {{T_{lk}}\vec{\boldsymbol\psi }_k^{nc}} } \right) \cdot \left( {\nabla \times \sum\limits_m {{T_{jm}}\vec{\boldsymbol\psi }_m^{nc}} } \right){\rm{d}}\Omega}. \end{eqnarray*} Now the summation over k and m can be taken off the sign of integration:   \begin{equation} {G_{lj}} = \frac{1}{{{\mu _0}}}\sum\limits_e {\sum\limits_k {\sum\limits_m {{T_{lk}}{T_{jm}}} } \int_{{{\Omega _e}}}{}\left( {\nabla \times \vec{\boldsymbol\psi }_k^{nc}} \right) \cdot \left( {\nabla \times \vec{\boldsymbol\psi }_m^{nc}} \right){\rm{d}}\Omega } . \end{equation} (63)Obviously, the integral in (63) defines the matrix $${{{\bf \hat{G}}}^e}$$ of the FE Ωe for basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$ which are the standard shape vector functions on Ωe. Thus, formula (63) defines the rule of assembling matrices $${{{\bf \hat{G}}}^e}$$ when the ‘conforming’ basis functions $$\vec{\boldsymbol\psi }_j^c$$ are used on the non-conforming meshes, and shape functions on FEs (and, consequently, their matrices) are defined in a standard way for the vector FEM. This rule is as follows. The numbers k and m of global basis functions from the basis $$\{ {\vec{\boldsymbol\psi }_i^{nc}} \}$$ are defined by the numbers of local basis functions, and the considered component of the local matrix is added with the factor TlkTjm only to those components Glj of the global matrix G, for which both components Tlk and Tjm are non-zero. The generation of matrices Mσ and Mε of FE equation system (56) is performed similarly, that is, with the use of standard shape vector functions, which correspond to the ‘non-conforming’ basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. Because the basis functions $$\vec{\boldsymbol\psi }_j^c$$ are the elements of the space Hcurl(Ω), all theoretical estimations of the convergence of FE solutions obtained with the use of these functions (Monk 2003) are true. If we perform the local mesh refinements in the needed subdomains (in the vicinity of the receivers and anomalies with significant influence), FE solutions obtained on such non-conforming meshes have essentially the same accuracy as FE solutions obtained on regular meshes, which correspond to these irregular non-conforming meshes. However, when such irregular meshes are used, the computational costs are lowered by an order or more (frequently, even more than 100 times) in comparison with the use of corresponding (i.e. with the cells refinement in the same subdomains) regular meshes. This will be shown below, in Section 5. 4.3 The grouping of logging tool positions and using the direct solver The application of the primary-secondary field approach with a separate calculation of the horizontally layered medium field and application of the unstructured optimized meshes with local refinements of the cells for calculating the secondary field allow us to significantly reduce the size of matrix of the FE system. In this case, the direct solvers become sufficiently effective. Because the solution of logging problems implies the computation of the large number of fields in the same medium for the different positions of transmitters and receivers, the application efficiency of the direct solvers can be highly increased. If the medium does not change and only the transmitter position changes, then the matrix of the FE system will remain the same for all these problems when the same mesh is used. In this case, only the right-hand side vectors of the FE system will be different (in accordance with the transmitter position). Thus, when the direct solvers for SLAE are used, the matrix factorization into triangular matrices can be performed only once for all transmitter positions, and when the field for each transmitter position is computed, one should solve only the systems with these triangular matrices (Da Silva et al.2012; Oldenburg et al.2013; Haber & Schwarzbach 2014). Note that in our code, we use Intel MKL PARDISO solver (Schenk & G̈artner 2004). However, simultaneous field calculation for all logging tool positions on the trajectory is inappropriate. First, the trajectories in the logging problems usually significantly exceed the size of the computational domain, in which the field for one transmitter position can be computed with a required accuracy. Second, in order to obtain the numerical solution with the required accuracy, it will suffice to refine the mesh only in a reasonably small vicinity of the tool position. Therefore, in order to achieve the required accuracy of the solution on the same mesh for several logging tool positions, it is required to refine the mesh in quite a large subdomain of the computational domain, and the number of nodes in the mesh (and, consequently, the SLAE size) can increase very much with the increase in the number of tool positions for which the same matrix of the FE system will be generated. As a result, the size of the FE system can increase so much that even the time of the solution of two systems with two triangular matrices (i.e. factors obtained due to the factorization) can exceed the total time of calculating the field for one tool position on the mesh with a small number of FEs (not to mention other computational costs: factorization of the big size matrix, recalculation of the primary field onto a fine mesh, and etc.). Thus, both the approach in which the field calculation is performed for each logging tool position on the individual small mesh and the approach in which the fields for all logging tool positions on the trajectory are calculated simultaneously (on the corresponding mesh with a great number of elements) are not optimal. Therefore, the most appropriate is the grouping of problems, which correspond to different tool positions on the trajectory, such that, on the one hand, there are more positions in one group, and, on the other hand, the size of the matrix of the FE system does not increase too much, because with the increase in the size of SLAE, the computational costs for matrix factorization (in order to construct the triangular factors) and other overhead costs grow very fast. In the case of increasing the number of positions in the group, the difference between the 3-D medium and the horizontally layered background medium taken for this group can also increase. In this case, in order to achieve the required accuracy, it is necessary to use a more fine mesh for calculating the secondary field that can lead to higher computational costs. We also note that the field calculation for groups can be performed in parallel that makes it possible to achieve the additional acceleration. In our code, the grouping of tool positions is performed automatically. In this case, the number of positions in a group is essentially defined by the number of positions getting into the sequent segment of the trajectory. The length of the trajectory segment (as well as the step sizes in the mesh) depends on the frequency. Thus, for the frequency of 100 kHz, the optimal length of the trajectory segment varies from 2 m (in the case, if the distance between adjacent positions is about 0.1 m) to 6 m (in the case, if the distance between adjacent positions is about 1–2 m). For other frequencies, the length of the trajectory segment is scaled by multiplying by the square root of the frequencies ratio. Thus, the number of positions in one group varies on the average from 3 to 20. 4.4 Automatic choice of the horizontally layered background medium for calculating the primary field The highest accuracy of the FE solution can be achieved without increasing the computational costs in the situation when the logging tool (i.e. its transmitter and receivers) is entirely located in such a subdomain of the 3-D computational domain the conductivity of which coincides with the conductivity of the background medium (i.e. the horizontally layered medium taken as a background for calculating the primary field). Moreover, the structures of geoelectrical models and logging tool motion trajectories, which are specific for the problems considered here (see e.g. Fig. 3), as a rule, make practically impossible the choice (pre-setting) of the single (for the entire trajectory) background medium, which would satisfy the requirements mentioned above. Therefore, we propose to perform the automatic choice of the background medium for each logging tool position (or group of positions) on the trajectory taking into account the location of the tool and the properties of the medium in close vicinity of it. For calculating the field for one group of logging tool positions, we propose the following algorithm of the automatic choice of the background medium. In the 3-D FE mesh, we select the subdomain ΩT-R, which in plane contains all transmitters and receivers of the considered group and some vicinity surrounding them (i.e. this subdomain is limited only in x- and y-directions, but in the z-direction, it coincides with the original 3-D computational domain). Using all FEs, which are inside this subdomain, the array ZT-R of z-coordinates of these FEs boundaries is generated. This array is ordered, and the numbers of FEs, the z-coordinates of which are located inside the interval ($$Z_i^{{\rm{T - R}}}$$, $$Z_{i + 1}^{{\rm{T - R}}}$$), are assigned to this interval. This information is formed as a list by one looping over all FEs from ΩT-R. Then, by means of looping over intervals the boundaries of which are defined by the array ZT-R, the horizontally layered (background) medium is formed by the following algorithm. For the current interval, we define the maximal conductivity in this interval by analysing all FEs associated with this interval and taken from the corresponding list. If this conductivity coincides with the maximal conductivity obtained for the previous interval, then we pass to the next interval. If we obtain a different value of maximal conductivity, then the interface between the current and previous intervals becomes the interface of the horizontally layered (background) medium to be formed. As conductivities of this horizontally layered (background) medium we take maximal conductivity values, which are obtained during the determination of the z-coordinates of layers interfaces. Because the size of the 3-D computational domain is smaller than the size of the computational domain for calculating the primary field, outside ΩT-R, the background horizontally layered medium is formed using the data describing the original problem. At the end of the algorithm, we check the thickness of the layers. If the layers with a very small thickness are presented in the formed medium, then these layers are joined with the adjacent ones (from two adjacent layers we choose the layer which is thicker), that is, in fact, a small thickness layer is excluded from the medium, and the layer adjacent to it is properly extended. We note that the 3-D mesh, which is used for the automatic choice of the background medium, can partially describe the geometry of the 3-D problem to be solved. Before automatically choosing, it is better not to include small size objects into this mesh (e.g. a borehole if its effect on the measured signal is rather small). We include these objects into the 3-D mesh after the layers of the background horizontally layered medium have been formed, and their interfaces have been included into the mesh. Fig. 9 (right) shows the sections of the background mediums, which are formed with the use of the algorithm considered above for three groups of tool positions, including the situation when the trajectory crosses the fault. Figure 9. View largeDownload slide The vertical sections of meshes (left) generated for different groups of adjacent tool positions for the case when the background medium (right) is chosen automatically. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the ith transmitter. Figure 9. View largeDownload slide The vertical sections of meshes (left) generated for different groups of adjacent tool positions for the case when the background medium (right) is chosen automatically. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the ith transmitter. Fig. 9 shows the FE meshes generated for the same tool positions on the trajectory as shown in Figs 5 and 6. But in Fig. 9, these meshes have been generated for the case when the background medium is chosen automatically for each group of the tool positions. The comparison of the numerical results for the cases when the background medium is pre-set as single and when the background medium is chosen automatically (for each group of the adjacent tool positions) will be presented below, in Section 6. 5 VALIDATION In this section, we present the results of the verification of the mathematical models and numerical schemes proposed in this paper and described in the previous sections. We will call the field calculation performed with the use of mathematical models (3), (20), (28)–(30) the 2-D approximation because, although the medium is horizontally layered (i.e. really 1-D), in order to calculate the electromagnetic field induced by the current coil, we use the solution of one or several axisymmetrical (i.e. 2-D) problems in the 2-D computational domain in the (r, z)-coordinates. We will call the field calculation performed with the use of mathematical model (1) the 3-D approximation even in the case when the 2.5-D problem is solved, because in order to obtain its solution, we find the FE solution in the 3-D computational domain. In order to verify the method proposed in this paper for calculating the primary field in the horizontally layered medium, we use the Dipole1D code (Key 2009). Because the Dipole1D code allows the calculation of the field induced only by the electrical dipole, then in order to obtain the field induced by the current frame (loop), we sum the fields of four dipoles making up the current frame. As a testing model, we use the geoelectrical model shown in Fig. 10. Figure 10. View largeDownload slide Testing geoelectrical model 1 of the horizontally layered medium and locations of the lines (Line1…Line5) along which the non-zero components of $${{\bf \vec{B}}}$$ and $${{\bf \vec{E}}}$$ are calculated with the use of the proposed approach and Dipole1D code (Key 2009). Figure 10. View largeDownload slide Testing geoelectrical model 1 of the horizontally layered medium and locations of the lines (Line1…Line5) along which the non-zero components of $${{\bf \vec{B}}}$$ and $${{\bf \vec{E}}}$$ are calculated with the use of the proposed approach and Dipole1D code (Key 2009). The source in the form of a square frame with the side of 0.1 m is fixed, and its centre is located at the point (0, 0, 13). The current is 100 A. The frequency is 100 kHz. In order to verify the field values in the receivers, we will output the non-zero components of the magnetic induction $${{{\bf \vec{B}}}^p}$$ along Line 1 and Line 2 (see Fig. 10). Besides, as we have described in the previous sections, when calculating the 3-D field, the right-hand side of the equations (see eq. 1) is defined by the electrical field induced by the source in the horizontally layered medium, which has been taken as a background. Therefore, we will also verify the accuracy of calculating the non-zero components of the electrical field $${{{\bf \vec{E}}}^p}$$. We will analyse the accuracy of calculating $${{{\bf \vec{E}}}^p}$$ along Line 3, Line 4, and Line 5 (see Fig. 10), and also along Line 6 and Line 7, which are parallel to Line 3 and have the same z-coordinate. The distance along the y-axis between these lines (Line 6 and Line 7) and Line 3 is 2 m and 4 m, respectively. The plots placed at the top of Fig. 11 show the real and imaginary x- and z-components of $${{{\bf \vec{B}}}^p}$$ along Line 1 and Line 2 for the case when the field is induced by the ‘vertical’ frame and ‘horizontal’ loop. The plots placed at the bottom of Fig. 11 show the relative (in per cent) deviations of $${{{\bf \vec{B}}}^p}$$, which are calculated using the approach proposed in this paper, from the values of $${{{\bf \vec{B}}}^p}$$, which are calculated with the use of the Dipole1D code (Key 2009). The results presented in Fig. 11 show that within the distance of 10 m, the relative deviations for all components do not exceed 0.5 per cent, and within the distance of 20 m, the deviations do not exceed 1–1.5 per cent. The error growth takes place in the subintervals of low values of the field. We do not consider the field values farther than 20 m, because this value (20 m) defines the maximal size of the 3-D computational domain. Figure 11. View largeDownload slide $$B_x^{{\rm{Re}}}$$ (dark blue curve), $$B_x^{{\mathop{\rm Im}\nolimits} }$$ (green curve), $$B_z^{{\rm{Re}}}$$ (red curve) and $$B_z^{{\mathop{\rm Im}\nolimits} }$$ (light blue curve) along Line 1 and Line 2 (see Fig. 10) (top) calculated using the method proposed in this paper and their deviations (bottom) from those that are obtained using the Dipole1D code (Key 2009) when the electromagnetic field is induced by the ‘vertical’ frame and ‘horizontal’ loop. Figure 11. View largeDownload slide $$B_x^{{\rm{Re}}}$$ (dark blue curve), $$B_x^{{\mathop{\rm Im}\nolimits} }$$ (green curve), $$B_z^{{\rm{Re}}}$$ (red curve) and $$B_z^{{\mathop{\rm Im}\nolimits} }$$ (light blue curve) along Line 1 and Line 2 (see Fig. 10) (top) calculated using the method proposed in this paper and their deviations (bottom) from those that are obtained using the Dipole1D code (Key 2009) when the electromagnetic field is induced by the ‘vertical’ frame and ‘horizontal’ loop. The plots at the top of Figs 12 and 13 show the non-zero components of $${{{\bf \vec{E}}}^p}$$ along Lines 3–7 for the case when the field induced by the ‘vertical’ frame (Fig. 12) and ‘horizontal’ loop (Fig. 13). The plots at the bottom of Figs 12 and 13 show the relative (in per cent) deviations of $${{{\bf \vec{E}}}^p}$$, which are calculated using the approach proposed in this paper, from the values of $${{{\bf \vec{E}}}^p}$$, which are calculated with the use of the Dipole1D code (Key 2009). The results presented in Figs 12 and 13 show that the error of the calculation does not essentially exceed 0.5 per cent. The error grows to 1–1.5 per cent only in the subintervals of low values of the field, which are not crucial for creating the 3-D field and, therefore, have practically no effect on its calculation accuracy. Figure 12. View largeDownload slide Non-zero components of $${{{\bf \vec{E}}}^p}$$ along lines 3–7 (see Fig. 10) calculated using the method proposed in this paper (top) and their deviations (bottom) from those that are obtained using the Dipole1D code (Key 2009) when the electromagnetic field is induced by the ‘vertical’ frame. Figure 12. View largeDownload slide Non-zero components of $${{{\bf \vec{E}}}^p}$$ along lines 3–7 (see Fig. 10) calculated using the method proposed in this paper (top) and their deviations (bottom) from those that are obtained using the Dipole1D code (Key 2009) when the electromagnetic field is induced by the ‘vertical’ frame. Figure 13. View largeDownload slide Non-zero components of $${{{\bf \vec{E}}}^p}$$ along lines 3–7 (see Fig. 10) calculated using the method proposed in this paper (top) and their deviations (bottom) from those that are obtained using the Dipole1D code (Key 2009) when the electromagnetic field is induced by the ‘horizontal’ loop. Figure 13. View largeDownload slide Non-zero components of $${{{\bf \vec{E}}}^p}$$ along lines 3–7 (see Fig. 10) calculated using the method proposed in this paper (top) and their deviations (bottom) from those that are obtained using the Dipole1D code (Key 2009) when the electromagnetic field is induced by the ‘horizontal’ loop. Thus, the obtained results verify the accuracy of calculating the electromagnetic field in the horizontally layered medium with the use of the method proposed in this paper. It is necessary to note that the number of points at which the output of $${{{\bf \vec{E}}}^p}$$ for the following solution of the 3-D problem is required is about a few tens of thousands. For example, for calculating the field at 30 000 points using the Dipole1D code, it is required about 90–130 s depending on the number of layers in the geoelectrical model. In this case, the computational time, which is required for calculating the field with the use of the method proposed in this paper, is about 10–12 s. Even if we assume that the computational time of the semi-analytical method can decrease by four times (we use four electrical dipoles for calculating the field induced by one source in the form of a frame or loop, because we have no software analogous to the Dipole1D code for the magnetic dipole), nevertheless, when calculating the primary field at 30 000 points, the method proposed in this paper makes it possible to reduce the computational time by more than two times. In addition, we note that for calculating the field in the group consisting of six positions of the tool, the number of points at which the primary field should be calculated can be about 100 000, and in this case, the computational advantage of the FE method proposed in this paper for calculating the primary field can reach an order. Therefore, the method of FE solution for the model proposed in this paper (Section 2.2) is much less expensive in comparison with the semi-analytical methods and, in addition, its accuracy is good enough. In order to verify the method of calculating the secondary field in the 3-D medium, we will compare the numerical results obtained for the horizontally layered model and curvilinear trajectory (Fig. 14) with the use of two approaches. In the first case, we will use the proposed method of the calculation in the horizontally layered medium (the 2-D approximation). In this case, the conductivity tensor is diagonal and has the same components (isotropic medium), which are defined as follows: $$\sigma _{ii}^p = 0.25$$ Sm m−1 if z < 490 m, $$\sigma _{ii}^p = 0.125$$ Sm m−1 if 490 m < z < 502 m, $$\sigma _{ii}^p = 0.025$$ Sm m−1 if 502 m < z < 507 m, $$\sigma _{ii}^p = 1$$ Sm m−1 if 507 m < z < 515 m, $$\sigma _{ii}^p = 0.333$$ Sm m−1 if z > 515 m. In the second case, we will use the method of 3-D field calculation (the 3-D approximation). In this case, the second, third, and fourth layers of the geoelectrical model considered above will be taken as 3-D geoelectrical inhomogeneities. The tensor σ is defined similarly to σp considered above, and σp in this problem is defined as: $$\sigma _{ii}^p = 0.25$$ Sm m−1 if z < 510 m, $$\sigma _{ii}^p = 0.333$$ Sm m−1 if z > 510 m. Figure 14. View largeDownload slide Testing geoelectrical model 2 of the horizontally layered medium. The transmitter and receivers are moved along the trajectory shown by solid black line. Figure 14. View largeDownload slide Testing geoelectrical model 2 of the horizontally layered medium. The transmitter and receivers are moved along the trajectory shown by solid black line. We will analyse the TZ-RZ array (T denotes a transmitter, R denotes a receiver) where Z is the direction of the logging tool axis, which coincides with its motion trajectory. The offsets of receivers from the transmitter are 6 m and 10 m, the frequency is 100 kHz. The distance between adjacent positions on the trajectory is 1 m. Figs 15 and 16 (left) show the amplitudes in two receivers (Fig. 15) and phase differences between two receivers offsets (Fig. 16) calculated with the use of 2-D approximation (corresponding curves are shown by dark blue colour) and 3-D approximation on two meshes (corresponding curves are shown by green colour for the coarse mesh and by red colour for the fine mesh). Fig. 15 (right) shows the relative deviations (in per cent) of the amplitudes obtained on the two meshes using the second approach (i.e. by solving the problem considered above using the 3-D approximation) from the solution obtained using the first approach (i.e. by solving the problem considered above using the 2-D approximation). The results presented in Fig. 15 show that the problem solution obtained with the use of 3-D approximation converges to the problem solution obtained with the use of 2-D approximation (the error of which does not exceed 0.5 per cent for offsets of 6 m and 10 m that has been shown above in Fig. 11). The error in the amplitudes obtained with the use of 3-D approximation does not exceed 2 per cent on the fine mesh. Fig. 16 (right) shows the absolute deviations (in degrees) of the phase differences obtained on the two meshes with the use of 3-D approximation from the phase differences obtained with the use of 2-D approximation. The results presented in Fig. 16 show that the error of the problem solution obtained with the use of 3-D approximation does not exceed 1° on the fine mesh. Figure 15. View largeDownload slide Amplitudes versus tool positions on the trajectory (left) and relative deviations (right) of the amplitudes obtained with the use of 3-D approximation on two meshes from the ones obtained with the use of 2-D approximation (dark blue curves). Green curves correspond to the coarse 3-D mesh, and red curves correspond to the fine 3-D mesh. Figure 15. View largeDownload slide Amplitudes versus tool positions on the trajectory (left) and relative deviations (right) of the amplitudes obtained with the use of 3-D approximation on two meshes from the ones obtained with the use of 2-D approximation (dark blue curves). Green curves correspond to the coarse 3-D mesh, and red curves correspond to the fine 3-D mesh. Figure 16. View largeDownload slide Phase differences between two receivers offsets (6 and 10 m) versus tool position on the trajectory (left) and deviations of the phase differences (right) obtained with the use of 3-D approximation on two meshes from the ones obtained with the use of 2-D approximation (dark blue curves). The green curve corresponds to the coarse 3-D mesh, and red curve corresponds to the fine 3-D mesh. Figure 16. View largeDownload slide Phase differences between two receivers offsets (6 and 10 m) versus tool position on the trajectory (left) and deviations of the phase differences (right) obtained with the use of 3-D approximation on two meshes from the ones obtained with the use of 2-D approximation (dark blue curves). The green curve corresponds to the coarse 3-D mesh, and red curve corresponds to the fine 3-D mesh. Below, in Section 6, we will show that if we use the procedure of automatic choice of the background medium, then the coarse mesh is sufficient for calculating the 3-D field with the error limited by 1 per cent. The verification of the developed numerical scheme is also performed by the comparison of numerical results obtained on the non-conforming irregular meshes and on the corresponding (with the same refinements in the vicinities of receivers and 3-D anomalies with the significant influence) regular meshes for the 3-D geoelectrical model presented in Fig. 3. The comparison shows that the deviations of the signals at the receivers do not exceed 0.7 per cent, and when the regular meshes are used, the calculation time grows by more than two orders. This is because the number of nodes in the regular meshes is approximately 20–30 times more than in the optimized irregular meshes corresponding to the regular meshes (i.e. with the cells refinements in the same subdomains). These irregular meshes have been presented in Sections 4.1 and 4.4. In order to verify the results of the electromagnetic field calculation in the anisotropic medium, we compare the solutions obtained for two geoelectrical models which are shown in Fig. 17. The first model (shown in Fig. 17a) is similar to the model shown in Fig. 14 except for the third layer with anisotropic conductivity: $$\sigma _{11}^p = 0.256$$ Sm m−1, $$\sigma _{22}^p = 0.025$$ Sm m−1 if 502 m < z < 507m. In the second model (shown in Fig. 17b), this layer has been replaced by a thin layered structure with the alternating conductivity values of 0.5\0.0128Sm m−1. We consider three variants of this layered structure: the thickness of layers is 0.5 m, 0.25 m, and 0.125 m. As before, we analyse the TZ-RZ array (where Z is the direction of the logging tool axis) for two frequencies of 100 and 200 kHz. The offsets of the receivers from the transmitter are 6 m and 10 m, respectively. Figure 17. View largeDownload slide (a) The testing model of the horizontally layered medium in which the third layer conductivity is anisotropic (lateral conductivity is 0.256 Sm m−1 and vertical conductivity is 0.025 Sm m−1) and (b) the model with the thin layered structure equivalent to the anisotropic layer (in the third layer, thin layers with alternate conductivity 0.5 Sm m−1 and 0.0128 Sm m−1 are set). Figure 17. View largeDownload slide (a) The testing model of the horizontally layered medium in which the third layer conductivity is anisotropic (lateral conductivity is 0.256 Sm m−1 and vertical conductivity is 0.025 Sm m−1) and (b) the model with the thin layered structure equivalent to the anisotropic layer (in the third layer, thin layers with alternate conductivity 0.5 Sm m−1 and 0.0128 Sm m−1 are set). Figs 18 and 19 (left) show the amplitudes and Figs 20 and 21 (left) show the phase differences computed for anisotropic (dark blue curves) and thin layered (green curves with labels) geoelectrical models in comparison with the curves computed for two isotropic geoelectrical models. The conductivity of the third layer in the first isotropic model (lilac curves) has been taken equal to 0.025 Sm m−1 (i.e. equal to the vertical conductivity of the third layer of the anisotropic model), and in the second model (brown curves) this conductivity has been taken equal to 0.256 Sm m−1 (i.e. equal to the lateral conductivity of the third layer of the anisotropic model). The numerical results presented in Figs 18–21 show that the deviation between the curves obtained for the anisotropic and isotropic media is significant, and the curves obtained for anisotropic and thin layered media are in good agreement. Figure 18. View largeDownload slide Amplitudes (left) obtained for the anisotropic medium (dark blue curve), thin layered structure (green curve with labels) and two isotropic media (brown and lilac curves), and relative deviations (right) of amplitudes obtained for thin layered models with layers’ thickness of 0.5 m (red curve), 0.25 m (light blue curve) and 0.125 m (green curve with labels) from those obtained for the anisotropic model. Frequency is 100 kHz. Figure 18. View largeDownload slide Amplitudes (left) obtained for the anisotropic medium (dark blue curve), thin layered structure (green curve with labels) and two isotropic media (brown and lilac curves), and relative deviations (right) of amplitudes obtained for thin layered models with layers’ thickness of 0.5 m (red curve), 0.25 m (light blue curve) and 0.125 m (green curve with labels) from those obtained for the anisotropic model. Frequency is 100 kHz. Figure 19. View largeDownload slide Amplitudes (left) obtained for the anisotropic medium (dark blue curve), thin layered structure (green curve with labels) and two isotropic media (brown and lilac curves), and relative deviations (right) of the amplitudes obtained for thin layered models with layers’ thickness of 0.5 m (red curve), 0.25 m (light blue curve) and 0.125 m (green curve with labels) from those obtained for the anisotropic model. Frequency is 200 kHz. Figure 19. View largeDownload slide Amplitudes (left) obtained for the anisotropic medium (dark blue curve), thin layered structure (green curve with labels) and two isotropic media (brown and lilac curves), and relative deviations (right) of the amplitudes obtained for thin layered models with layers’ thickness of 0.5 m (red curve), 0.25 m (light blue curve) and 0.125 m (green curve with labels) from those obtained for the anisotropic model. Frequency is 200 kHz. Figure 20. View largeDownload slide Phase differences (left) obtained for the anisotropic medium (dark blue curve), thin layered structure (green curve with labels) and two isotropic media (brown and lilac curves), and deviations (right) of the phase differences obtained for thin layered models with layers’ thickness of 0.5 m (red curve), 0.25 m (light blue curve) and 0.125 m (green curve with labels) from those obtained for the anisotropic model. Frequency is 100 kHz. Figure 20. View largeDownload slide Phase differences (left) obtained for the anisotropic medium (dark blue curve), thin layered structure (green curve with labels) and two isotropic media (brown and lilac curves), and deviations (right) of the phase differences obtained for thin layered models with layers’ thickness of 0.5 m (red curve), 0.25 m (light blue curve) and 0.125 m (green curve with labels) from those obtained for the anisotropic model. Frequency is 100 kHz. Figure 21. View largeDownload slide Phase difference (left) curves obtained for the anisotropic medium (dark blue curve), thin layered structure (green curve with labels) and two isotropic media (brown and lilac curves), and deviations (right) of the phase differences obtained for thin layered models with layers’ thickness of 0.5 m (red curve), 0.25 m (light blue curve) and 0.125 m (green curve with labels) from those obtained for the anisotropic model. Frequency is 200 kHz. Figure 21. View largeDownload slide Phase difference (left) curves obtained for the anisotropic medium (dark blue curve), thin layered structure (green curve with labels) and two isotropic media (brown and lilac curves), and deviations (right) of the phase differences obtained for thin layered models with layers’ thickness of 0.5 m (red curve), 0.25 m (light blue curve) and 0.125 m (green curve with labels) from those obtained for the anisotropic model. Frequency is 200 kHz. Figs 18–21 (right) show the deviations of the curves obtained for three variants of the model in which the anisotropic layer represented by the thin layered structure with the thickness of 0.5 m (red curve), 0.25 m (light blue curve), and 0.125 m (green curve with labels) from the curves obtained for the geoelectrical model with anisotropic conductivity of this layer. The numerical results presented in Figs 18–21 (right) show that when the layers thickness decreases (i.e. the number of layers increases), then the amplitudes and phase differences obtained for the thin layered medium approach the amplitudes and phase differences obtained for the anisotropic model. These results, on the one hand, confirm the correctness of the implementations performed taking into account the anisotropic conductivity, and, on the other hand, make it possible to define the layers thickness starting with which the thin layered structure can be replaced by a layer with anisotropic conductivity. Note that for the tool parameters (and the frequencies used in them) which have been considered above, the effect of displacement currents is insignificant. Therefore, in order to verify that these currents are correctly accounted, we use another tool with a much higher frequency. We perform the calculations for two variants of the medium shown in Fig. 22(a). In the first variant, we set the relative dielectric permeability $$\tilde{\varepsilon }$$ of the low-conductive layer equal to 10 (see Fig. 22a, model 1), and in the second variant, we set its value equal to 1 (see Fig. 22a, model 2). As before, we analyse the TZ-RZ array. The frequency is equal to 14 MHz, and the offset of receiver from the transmitter is 0.5 m. The tool moves along the trajectory perpendicular to the layers of the medium. We calculate the field using two approaches. In the first approach, we use the 3-D approximation, that is, the mathematical model (1) in which σ = σp for the geoelectrical model under consideration, and, in this case, the secondary field is determined by the displacement current effect. In the second approach, we use the 2-D approximation, that is, the equation   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta E_\varphi ^{} + \frac{1}{{{\mu _0}{r^2}}}E_\varphi ^{} + i\omega \sigma _r^{p,rz}E_\varphi ^{} - {\omega ^2}\varepsilon E_\varphi ^{} = - i\omega I{\delta ^{{J_\varphi }}}. \end{equation} (64) Figure 22. View largeDownload slide Verification of numerical results obtained taking into account the displacement current effect for the frequency equal to 14 MHz and offset equal to 0.5 m: (a) geoelectrical model; (b) amplitudes calculated for model 1 without regard for the displacement current effect (green curve) and with regard for the displacement current effect when 2-D (dark blue curve) and 3-D (red curve) approximations are used; (c) amplitudes calculated for model 1 (dark blue curve) and for model 2 (red curve) in comparison with the ones calculated without regard for the displacement current effect (green curve). Figure 22. View largeDownload slide Verification of numerical results obtained taking into account the displacement current effect for the frequency equal to 14 MHz and offset equal to 0.5 m: (a) geoelectrical model; (b) amplitudes calculated for model 1 without regard for the displacement current effect (green curve) and with regard for the displacement current effect when 2-D (dark blue curve) and 3-D (red curve) approximations are used; (c) amplitudes calculated for model 1 (dark blue curve) and for model 2 (red curve) in comparison with the ones calculated without regard for the displacement current effect (green curve). We can use this equation, because we consider the true-vertical trajectory and the TZ-RZ array. Fig. 22(b) shows the amplitudes calculated for geoelectrical model 1 with the use of 2-D (dark blue curve) and 3-D (red curve) approximation. In this figure, for comparison, we show the amplitude, which is calculated without regard for the displacement current effect (green curve). One can see that the numerical results obtained with the use of 2-D and 3-D approximation are in good agreement (that confirms the correctness of implementations of 2-D and 3-D approximation and sufficient accuracy of 3-D calculations). Moreover, we can see the degree of displacement current effect for this logging tool and the geoelectrical model. Fig. 22(c) presents the comparison of amplitudes calculated for model 1 (dark blue curve) and model 2 (red curve), which differ from each other only in the values of $$\tilde{\varepsilon }$$ in the low-conductive layer with the conductivity equal to 0.025 Sm m−1. The results presented in Fig. 22 show that on the trajectory segment within this layer, the amplitude values obtained for model 2 (red curve) approach the values calculated without regard for the displacement current (green curve). In order to additionally verify our code, we compare the numerical results, which are obtained by the approaches and code proposed in this paper, with the results presented in the paper (Hou et al.2006), which are obtained with the use of FD and IE methods for the 3-D geoelectrical model shown in Fig. 23. This model represents the seven-layered medium (some layers of which are anisotropic) and the borehole with the invasion zone. The diameter of the borehole is 21.59 cm, its conductivity is 2.5 Sm m−1, and the dip angle is 60°. The invasion zone is in the isotropic layers, and its diameter is 30.48 cm. Figure 23. View largeDownload slide The geoelectrical model with anisotropic layers and the borehole with the invasion zone. The diameter of the borehole is 21.59 cm, its conductivity is 2.5 Sm m−1, and the dip angle is 60°. The invasion zone is in isotropic layers, its diameter is 30.48 cm and conductivity is 0.2 Sm m−1 in the second layer (shown by red colour), 0.5 Sm m−1 in the fourth layer (shown by yellow colour), and 2 Sm m−1 in the six layer (shown by blue colour). Figure 23. View largeDownload slide The geoelectrical model with anisotropic layers and the borehole with the invasion zone. The diameter of the borehole is 21.59 cm, its conductivity is 2.5 Sm m−1, and the dip angle is 60°. The invasion zone is in isotropic layers, its diameter is 30.48 cm and conductivity is 0.2 Sm m−1 in the second layer (shown by red colour), 0.5 Sm m−1 in the fourth layer (shown by yellow colour), and 2 Sm m−1 in the six layer (shown by blue colour). Fig. 24 shows the imaginary components of the magnetic field Hxx for the frequencies of 100 kHz and 50 kHz when the transmitters and receivers coils axis are orthogonal to the borehole axis. The offset of the receiver from the transmitter is 1.8 m. The curves, which are obtained by our code, are shown by solid dark blue lines, and the results, which were obtained by FD and IE methods and presented in (Hou et al.2006), are shown by green and red points, respectively. Besides, Fig. 24 shows the curves calculated for the 1-D model (i.e. without the borehole and invasion zone). Figure 24. View largeDownload slide The imaginary components of the magnetic field Hxx for the frequencies of (a) 100 kHz and (b) 50 kHz when the transmitters’ and receivers’ coil axes are orthogonal to the borehole axis. The offset of the receiver from transmitter is 1.8 m. The curves, which are obtained by our code, are shown by solid dark blue lines, and the results, which were obtained by FD and IE methods and presented in Hou et al. (2006), are shown by green and red points, respectively. The curves calculated for the 1-D model (i.e. without the borehole and invasion zone) are shown by brick-red lines. Figure 24. View largeDownload slide The imaginary components of the magnetic field Hxx for the frequencies of (a) 100 kHz and (b) 50 kHz when the transmitters’ and receivers’ coil axes are orthogonal to the borehole axis. The offset of the receiver from transmitter is 1.8 m. The curves, which are obtained by our code, are shown by solid dark blue lines, and the results, which were obtained by FD and IE methods and presented in Hou et al. (2006), are shown by green and red points, respectively. The curves calculated for the 1-D model (i.e. without the borehole and invasion zone) are shown by brick-red lines. The presented results show that, on the whole, the values of Hxx calculated by our code lie between the values obtained by means of FD and IE methods (i.e. quite close to both results). Note that in the zone of intersection of the low-conductivity layer (z is from 0 to 1.5 m), our data (values of Hxx) are in good agreement with FD-results, and in the zone of intersection of the high-conductivity layer (z is from 7.875 m to 8.75 m), our data are in good agreement with IE-results. The modelling is performed for 118 tool positions (the step between the adjacent positions is 0.2 m) for two frequencies, that is, 236 problems are solved. The computational time for all positions and two frequencies is 12 min 15 s, that is, 3.11 s per point, if one core of an Intel Core i7–3770 K CPU 3.5 GHz is used. The memory required does not exceed 1 GB. When two, three, and four cores are used, the computational time for all positions and two frequencies is 6 min 10 s, 4 min 8 s, and 3 min 25 s (i.e. 1.57 s, 1.05 s, and 0.87 s per point), respectively. The results given in (Hou et al.2006) were obtained with the use of Intel Pentium 4 CPU 3.2 GHz (its calculation capacity is about five times lower than that of the computer used for our calculations). The computational time presented in (Hou et al.2006) for 50 positions is 1.9 hr in the sequential mode that corresponds to 2.28 min per point. Thus, in the sequential mode, the computational time per point required for our code, when the Intel Core i7–3770 K CPU 3.5 GHz is used, is about forty five times less than that presented in (Hou et al.2006) (i.e. when the calculations were performed by the FD-approach proposed in (Hou et al.2006) on an Intel Pentium 4 CPU 3.2 GHz). If we recalculate the computational time which would be required for calculating on the computer with the same capacity, then our code performs the calculation per point about ten times faster. The FE mesh used for calculating the 3-D (secondary) field for a group of tool positions (approximately 15 positions in one group) contains about 45 000 elements (that corresponds to about 110 000 complex-valued degrees of freedom). The 2-D meshes used for calculating the primary field contain about 20 000 nodes. 6 NUMERICAL RESULTS In this section, we present the numerical results obtained for the geoelectrical model with non-parallel layers and the fault, which is shown in Fig. 3. We consider the anisotropic model (the conductivity values of structural parts of this model are shown in Fig. 3) and two isotropic models. In one of these models, the conductivity of the third layer is set equal to the vertical conductivity (σ = 0.025 Sm m−1, further, we will call it ‘isotropic model 1’), and in the second model, the conductivity of the third layer is set equal to the lateral conductivity (σ = 0.256 Sm m−1, further, we will call it ‘isotropic model 2’). First, we present the comparison of numerical results obtained when the background medium is fixed (pre-set) for the entire trajectory, and when the background medium is automatically chosen during the calculation process (see Section 4.4). Figs 25 and 26 (left) show the amplitudes in two receivers (with offsets of 6 m and 10 m from the transmitter) and phase differences between their offsets calculated for isotropic model 1 along the trajectory shown in Fig. 3. The calculations have been performed on two meshes with different numbers of nodes and for two variants: when the background medium is fixed for the entire trajectory (variant 1) and when the background medium is automatically chosen for each group of tool positions (variant 2). In variant 1, the fixed background medium contains two layers with the conductivities of 0.25 Sm m−1 and 0.33 Sm m−1; the Z-coordinate of the interface of these layers is 510 m. The fine mesh has been obtained from the coarse mesh by double decreasing the step in three directions in the base mesh. The length of the trajectory segment for one group has been taken equal to 6 m. Because the step along the trajectory has been taken equal to 1 m, each group contains 6 positions of the logging tool. The coarse mesh for one group contains, on average, about 25 000 elements (that corresponds to about 70 000 complex-valued degrees of freedom), and the fine mesh contains about 55 000 elements (that corresponds to about 140 000 complex-valued degrees of freedom). Figs 25 and 26 (right) show the deviations of the amplitudes and phase differences calculated on the coarse meshes for variant 1 (red curves) and variant 2 (green curves), and on the fine mesh for variant 1 (dark blue curves) from the amplitudes and phase differences calculated on the fine mesh for variant 2. Figure 25. View largeDownload slide Amplitudes in two receivers obtained for isotropic model 1 (left) and relative deviations (right) of the amplitudes calculated for the fixed background medium on the coarse (red curve) and fine (dark blue curve) meshes and with the automatic choice of the background medium on the coarse mesh (green curve) from the amplitudes calculated with the automatic choice of the background medium on the fine mesh. Frequency is 100 kHz. Figure 25. View largeDownload slide Amplitudes in two receivers obtained for isotropic model 1 (left) and relative deviations (right) of the amplitudes calculated for the fixed background medium on the coarse (red curve) and fine (dark blue curve) meshes and with the automatic choice of the background medium on the coarse mesh (green curve) from the amplitudes calculated with the automatic choice of the background medium on the fine mesh. Frequency is 100 kHz. Figure 26. View largeDownload slide Phase differences in two receivers obtained for isotropic model 1 (left) and the deviations (right) of phase differences calculated for the fixed background medium on the coarse (red curve) and fine (dark blue curve) meshes and with the automatic choice of the background medium on the coarse mesh (green curve) from the phase differences calculated with the automatic choice of the background medium on the fine mesh. Frequency is 100 kHz. Figure 26. View largeDownload slide Phase differences in two receivers obtained for isotropic model 1 (left) and the deviations (right) of phase differences calculated for the fixed background medium on the coarse (red curve) and fine (dark blue curve) meshes and with the automatic choice of the background medium on the coarse mesh (green curve) from the phase differences calculated with the automatic choice of the background medium on the fine mesh. Frequency is 100 kHz. The results presented in Figs 25 and 26 show that in variant 2, the sufficiently accurate values of the electromagnetic field characteristics can be obtained even on the coarse mesh (the relative deviations, on average, do not exceed 1–1.5 per cent in the amplitudes and 1–1.5 degrees in the phase difference). In variant 1, the field values in the receivers have been calculated with the error reaching 6 per cent in the amplitudes and 3 degrees in the phase difference on the coarse mesh and with the error reaching 1–1.5 per cent in the amplitudes and 1–1.5 degrees in the phase difference on the fine mesh. In addition, we note that the meshes used in variant 1 are more refined than the meshes used in variant 2. The computational time for 152 logging tool positions and two frequencies (i.e. 304 3-D problems are solved) on the coarse mesh in variant 2 is 15 min 30 s on one computational core (i.e. in the sequential mode) of Intel Core i7–3770 K 3.5 GHz. When we use parallelization on 2, 3 and 4 cores, the computational time is 7 min 40 s, 5 min 5 s and 4 min, respectively. Each group contains about six tool positions on the trajectory. In variant 1 (i.e. without automatically choosing the background medium), the computational time required for achieving the equivalent accuracy of the solution is 12 min when we use four computational cores. Thus, the proposed algorithm of automatically choosing the background medium makes it possible to reduce the computational costs by almost 3 times. This advantage will increase with the growth of the contrast of the conductivities of the structural parts of the geoelectrical model. It should be noted again that we managed to provide such relatively low computational costs for solving quite a large set of problems not only due to the use of the non-conforming mesh for the FE approximation of the 3-D problem but also due to the use of the FE approach proposed in this paper (Section 2.2) for calculating the primary field in the 1-D background medium. For the solution of 304 3-D problems, it has been required to calculate the values of the primary field at approximately one million edges of FEs which lie inside the 3-D anomalies (in the coarse mesh). This means that, when the Dipole1D code (Key 2009) is used with parallelization on 4 cores, the computational costs required for calculating the primary field will be about 15 min. These costs are many times higher than the costs required for solving a whole set of 3-D problems with the use of the approach proposed in this paper including the field calculation in the background medium. Compare the amplitudes and phase differences along the trajectory for the anisotropic model (see Fig. 3) and two isotropic models considered above. The amplitudes and phase differences in two receivers (6 and 10 m) for the TZ-RZ array and for two frequencies (100 and 200 kHz) are shown in Figs 27 and 28. The places where the trajectory intersects the borders between the layers with different conductivities are shown by vertical grey dashed lines. Figure 27. View largeDownload slide Amplitudes obtained for TZ-RZ array for the anisotropic model (dark blue curves), isotropic model 1 (lilac curves) and isotropic model 2 (brown curves) for two frequencies: 100 kHz (left) and 200 kHz (right). The places where the trajectory intersects the borders between the layers with different conductivities are shown by vertical grey dashed lines. Figure 27. View largeDownload slide Amplitudes obtained for TZ-RZ array for the anisotropic model (dark blue curves), isotropic model 1 (lilac curves) and isotropic model 2 (brown curves) for two frequencies: 100 kHz (left) and 200 kHz (right). The places where the trajectory intersects the borders between the layers with different conductivities are shown by vertical grey dashed lines. Figure 28. View largeDownload slide Phase differences obtained for TZ-RZ array for the anisotropic model (dark blue curves), isotropic model 1 (lilac curves) and isotropic model 2 (brown curves) for two frequencies: 100 kHz (left) and 200 kHz (right). The places where the trajectory intersects the borders between the layers with different conductivities are shown by vertical grey dashed lines. Figure 28. View largeDownload slide Phase differences obtained for TZ-RZ array for the anisotropic model (dark blue curves), isotropic model 1 (lilac curves) and isotropic model 2 (brown curves) for two frequencies: 100 kHz (left) and 200 kHz (right). The places where the trajectory intersects the borders between the layers with different conductivities are shown by vertical grey dashed lines. The numerical results presented in Figs 27 and 28 show that the curves obtained for the anisotropic model significantly differ from the curves obtained for isotropic model 2 (in which the conductivity of the third layer has been taken equal to the lateral conductivity of the third layer in the anisotropic model), but they are similar to the curves which are obtained for isotropic model 1 (in which the conductivity of the third layer has been taken equal to the vertical conductivity of the third layer in the anisotropic model). Fig. 29 shows the amplitudes in the first receiver (with the offset equal to 6 m) for the TX-RZ array, where X is the direction perpendicular to the tool axis and lies in the plane of the geoelectrical model section shown in Fig. 3. We can see that for this array (in comparison with the TZ-RZ array) the curves calculated for the anisotropic model are very similar to the curves calculated for isotropic model 2 and differ from the curves calculated for isotropic model 1. Figure 29. View largeDownload slide Amplitudes obtained for TX-RZ array for the anisotropic model (dark blue curves), isotropic model 1 (lilac curves) and isotropic model 2 (brown curves) for two frequencies: 100 kHz (left) and 200 kHz (right). The places where the trajectory intersects the borders between the layers with different conductivities are shown by vertical grey dashed lines. Figure 29. View largeDownload slide Amplitudes obtained for TX-RZ array for the anisotropic model (dark blue curves), isotropic model 1 (lilac curves) and isotropic model 2 (brown curves) for two frequencies: 100 kHz (left) and 200 kHz (right). The places where the trajectory intersects the borders between the layers with different conductivities are shown by vertical grey dashed lines. Thus, joint analysis of the amplitudes and phase differences for TZ-RZ and TX-RZ arrays makes it possible to define that the crossed layer has the anisotropic conductivity, and this property can be used when developing the inversion procedures. 7 CONCLUSIONS The approach proposed in this paper for solving multisource induction logging problems in multidimensional media makes it possible to significantly reduce the computational costs together with providing high accuracy of numerical solutions when the logging tool moves in the inclined and horizontal wells, and when the inclination angle of the well changes. The computational effectiveness of the approach proposed is ensured by three major factors. First, we use the primary-secondary field approach in which, for calculating the primary field, the method based on the FE solution to several axisymmetrical problems is proposed. When the 3-D problem is solved by FEM, the primary-secondary field approach by itself gives the possibility of reducing the computational costs by an order or more, especially, if the background medium is similar to the 3-D medium in which the field should be calculated. But this approach requires the calculation of the primary field at a great number of points, and the FE method proposed for calculating the primary field makes the primary-secondary field approach even more computationally effective. The semi-analytical methods developed to date are very effective for calculating the 1-D field at a small number of points. However, when it is necessary to calculate the field at a great number of points, their computational costs can become very high. For the examples considered in this paper, we show that the computational costs of calculating the total 3-D fields for one trajectory of the logging tool motion with the use of primary-secondary field approach and the method proposed for calculating the primary field are several times lower than the computational costs required for calculating only the primary field (at the required number of points) with the use of the Dipole1D code (Key 2009) that implements the semi-analytical method. Second, for calculating the secondary field, the irregular non-conforming meshes (with hanging nodes) are used that allows local refinements of the mesh cells only in the required subdomains where it is necessary to approximate the boundaries of 2-D or 3-D anomalies in more detail or calculate the field with higher accuracy. This makes it possible to reduce the number of degrees of freedom (the unknowns) in the problem under consideration by ten or more times in comparison with the corresponding regular meshes. The generation of the irregular non-conforming meshes is performed fully automatically taking into account not only the boundaries of 2-D or 3-D anomalies but also the chosen background medium (the most similar to the 3-D medium where the processed group of the logging tool positions is located). For such non-conforming meshes, we have developed the method of constructing the ‘conforming’ basis functions from the ‘non-conforming’ ones and the algorithm of generating the global matrices and right-hand side vectors of the FE equation system for this ‘conforming’ basis by assembling the matrices and vectors of FEs with standard local basis vector functions, which correspond to non-conforming global basis functions. In general, the application of the irregular non-conforming meshes in comparison with that of the corresponding regular meshes enables us to reduce the computational time required for the 3-D problem by hundred or more times without any loss of accuracy. Third, in order to solve the induction logging problem with quite a long trajectory and a great number of the logging tool positions, the method of grouping these positions has been proposed. It allows significantly increasing the computational effectiveness, when the direct solver for the FE system is used, due to the compromise between the decrease in the number of the FE matrix factorizations with the increasing number of positions in one group and a sharp decrease in the FE system size for a small group of positions if the automatically chosen (for this group) background medium is quite similar to the modelled 3-D medium (i.e. the primary field does not differ too much from the total field). This optimized grouping of the logging tool positions allows the additional reduction in the computational costs, at least, by several times in comparison with the approaches which do not use grouping (even if the iterative solver is used) or, on the contrary, use the grouping of a great number of the logging tool positions. The developed approaches have been verified in the conductivity range of 0.01 to 10 Sm m−1 and in the frequency range of 10 kHz to 14 MHz. The algorithms proposed for mesh generation provide the required solution accuracy, and the number of mesh cells does not exceed 70 000. This allows determining the fixed memory requirements. Thus, when the calculation is performed on one computational core, the required memory does not exceed 2 GB. The application of the direct solver and algorithms developed for mesh generation provide a small dependence of computational time on the conductivity contrasts and computational domain structure, and the algorithm of automatic choice of the background horizontally layered medium provides the required solution accuracy if the conductivity contrasts are within the range of 0.01÷10 Sm m−1 when the computational domain structure is quite complex. In general, when modelling the logging tool motion along an arbitrary trajectory, the approach proposed in this paper makes it possible to solve several hundred 3-D problems within several minutes using the PC with Quad core Intel i7 3.5 GHz processor. The developed algorithms for generating the FE SLAE on the non-conforming meshes can be useful for solving the marine and aero electromagnetic problems including the problems in time-domain. However, the principles of automatic generation of these meshes in marine and aero electromagnetic problems can be different, and therefore, the code needs to be modified. Acknowledgements The study was performed under the applied research with the financial support from the Ministry of Education and Science of the Russian Federation (project unique identifier: RFMEFI57716×0216, project No. 14.577.21.0216). REFERENCES Ansari S., Farquharson C.G., 2014. 3D finite-element forward modeling of electromagnetic data using vector and scalar potentials and unstructured grids, Geophysics , 79, E149– E165. https://doi.org/10.1190/geo2013-0172.1 Google Scholar CrossRef Search ADS   Badea E.A., Everett M.E., Newman G.A., Biro O., 2001. Finite-element analysis of controlled-source electromagnetic induction using Coulomb-gauged potentials, Geophysics , 66, 786– 799. https://doi.org/10.1190/1.1444968 Google Scholar CrossRef Search ADS   Börner R.U., Ernst O.G., Güttel S., 2015. Three-dimensional transient electromagnetic modelling using Rational Krylov methods, Geophys. J. Int. , 202, 2025– 2043. https://doi.org/10.1093/gji/ggv224 Google Scholar CrossRef Search ADS   Chung Y., Son J.S., Lee T.J., Kim H.J., Shin C., 2014. Three-dimensional modelling of controlled-source electromagnetic surveys using an edge finite-element method with a direct solver, Geophys. Prospect. , 62, 1468– 1483. https://doi.org/10.1111/1365-2478.12132 Google Scholar CrossRef Search ADS   Da Silva N.V., Morgan J.V., MacGregor L., Warner M., 2012. A finite element multifrontal method for 3D CSEM modeling in the frequency domain, Geophysics , 77, E101– E115. https://doi.org/10.1190/geo2010-0398.1 Google Scholar CrossRef Search ADS   Fang S., Gao G.-Z., Torres-Verdín C., 2003. Efficient 3-D electromagnetic modeling in the presence of anisotropic conductive media using integral equations, in Three-Dimensional Electromagnetics Workshop No. 3 (3DEM-3), Extended Abstracts , pp. 1– 8. Google Scholar CrossRef Search ADS   Farquharson C.G., Miensopust M.P., 2011. Three-dimensional finite-element modelling of magnetotelluric data with a divergence correction, J. Appl. Geophys. , 75, 699– 710. https://doi.org/10.1016/j.jappgeo.2011.09.025 Google Scholar CrossRef Search ADS   Grayver A.V., 2015. Parallel three-dimensional magnetotelluric inversion using adaptive finite-element method. Part I: Theory and synthetic study, Geophys. J. Int. , 202, 584– 603. https://doi.org/10.1093/gji/ggv165 Google Scholar CrossRef Search ADS   Grayver A.V., Bürg M., 2014. Robust and scalable 3-D geo-electromagnetic modelling approach using the finite element method, Geophys. J. Int. , 198, 110– 125. https://doi.org/10.1093/gji/ggu119 Google Scholar CrossRef Search ADS   Grayver A.V., Kolev T.V., 2015. Large-scale 3D geoelectromagnetic modeling using parallel adaptive high-order finite element method, Geophysics , 80, E277– E291. https://doi.org/10.1190/geo2015-0013.1 Google Scholar CrossRef Search ADS   Grayver A.V., Streich R., Ritter O., 2013. Three-dimensional parallel distributed inversion of CSEM data using a direct forward solver, Geophys. J. Int. , 193, 1432– 1446. https://doi.org/10.1093/gji/ggt055 Google Scholar CrossRef Search ADS   Haber E., Ruthotto L., 2014. A multiscale finite volume method for Maxwell's equations at low frequencies, Geophys. J. Int. , 199, 1268– 1277. https://doi.org/10.1093/gji/ggu268 Google Scholar CrossRef Search ADS   Haber E., Schwarzbach C., 2014. Parallel inversion of large-scale airborne time-domain electromagnetic data with multiple OcTree meshes, Inverse Probl. , 30( 5), 055011. https://doi.org/10.1088/0266-5611/30/5/055011 Google Scholar CrossRef Search ADS   Hennessy L., Macnae J., 2010. Inversion of concentric loop electromagnetic data by transformation to an equivalent potential field response, Explor. Geophys. , 41, 240– 249. https://doi.org/10.1071/EG10022 Google Scholar CrossRef Search ADS   Hou J., Mallan R.K., Torres-Verdín C., 2006. Finite-difference simulation of borehole EM measurements in 3D anisotropic media using coupled scalar-vector potentials, Geophysics , 71, G225– G233. https://doi.org/10.1190/1.2245467 Google Scholar CrossRef Search ADS   Hunkeler P.A. et al.  , 2016. Improved 1D inversions for sea ice thickness and conductivity from electromagnetic induction data: Inclusion of nonlinearities caused by passive bucking, Geophysics , 81, WA45– WA58. https://doi.org/10.1190/geo2015-0130.1 Google Scholar CrossRef Search ADS   Jahandari H., Farquharson C.G., 2014. A finite-volume solution to the geophysical electromagnetic forward problem using unstructured grids, Geophysics , 79, E287– E302. https://doi.org/10.1190/geo2013-0312.1 Google Scholar CrossRef Search ADS   Jahandari H., Farquharson C.G., 2015. Finite-volume modelling of geophysical electromagnetic data on unstructured grids using potentials, Geophys. J. Int. , 202, 1859– 1876. https://doi.org/10.1093/gji/ggv257 Google Scholar CrossRef Search ADS   Jaysaval P., Shantsev D.V., de la Kethulle de Ryhove S., 2015. Efficient 3-D controlled-source electromagnetic modelling using an exponential finite-difference method, Geophys. J. Int. , 203, 1541– 1574. https://doi.org/10.1093/gji/ggv377 Google Scholar CrossRef Search ADS   Kamm J., Pedersen L.B., 2014. Inversion of airborne tensor VLF data using integral equations, Geophys. J. Int. , 198, 775– 794. https://doi.org/10.1093/gji/ggu161 Google Scholar CrossRef Search ADS   Key K., 2009. 1D inversion of multicomponent, multifrequency marine CSEM data: Methodology and synthetic studies for resolving thin resistive layers, Geophysics , 74, F9– F20. https://doi.org/10.1190/1.3058434 Google Scholar CrossRef Search ADS   Kordy M., Wannamaker P., Maris V., Cherkaev E., Hill G., 2016a. 3-D magnetotelluric inversion including topography using deformed hexahedral edge finite elements and direct solvers parallelized on SMP computers - Part I: Forward problem and parameter Jacobians, Geophys. J. Int. , 204, 74– 93. https://doi.org/10.1093/gji/ggv410 Google Scholar CrossRef Search ADS   Kordy M., Wannamaker P., Maris V., Cherkaev E., Hill G., 2016b. 3-dimensional magnetotelluric inversion including topography using deformed hexahedral edge finite elements and direct solvers parallelized on symmetric multiprocessor computers - Part II: Direct data-space inverse solution, Geophys. J. Int. , 204, 94– 110. https://doi.org/10.1093/gji/ggv411 Google Scholar CrossRef Search ADS   Lelièvre P.G., Farquharson C.G., 2013. Gradient and smoothness regularization operators for geophysical inversion on unstructured meshes, Geophys. J. Int. , 195, 330– 341. https://doi.org/10.1093/gji/ggt255 Google Scholar CrossRef Search ADS   Ley-Cooper A.Y., Viezzoli A., Guillemoteau J., Vignoli G., Macnae J., Cox L., Munday T., 2015. Airborne electromagnetic modelling options and their consequences in target definition, Explor. Geophys. , 46, 74– 84. https://doi.org/10.1071/EG14045 Google Scholar CrossRef Search ADS   Macnae J., 2015. 3D-spectral CDIs: a fast alternative to 3D inversion?, Explor. Geophys. , 46, 12– 18. https://doi.org/10.1071/EG14036 Google Scholar CrossRef Search ADS   Mogilatov V., Goldman M., Persova M., Soloveichik Y., Koshkina Y., Trubacheva O., Zlobinskiy A., 2016. Application of the marine circular electric dipole method in high latitude Arctic regions using drifting ice floes, J. Appl. Geophys. , 135, 17– 31. https://doi.org/10.1016/j.jappgeo.2016.08.007 Google Scholar CrossRef Search ADS   Monk P., 2003. Finite Element Methods for Maxwell's Equations , Oxford University Press. Google Scholar CrossRef Search ADS   Mukherjee S., Everett M.E., 2011. 3D controlled-source electromagnetic edge-based finite element modeling of conductive and permeable heterogeneities, Geophysics , 76, F215–F226. https://doi.org/10.1190/1.3571045 Newman G.A., Alumbaugh D.L., 1995. Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences, Geophys. Prospect. , 43, 1021– 1042. https://doi.org/10.1111/j.1365-2478.1995.tb00294.x Google Scholar CrossRef Search ADS   Nie X.C., Yuan N., Liu R., 2013. A fast integral equation solver for 3D induction well logging in formations with large conductivity contrasts, Geophys. Prospect. , 61, 645– 657. https://doi.org/10.1111/j.1365-2478.2012.01070.x Google Scholar CrossRef Search ADS   Oldenburg D.W., Haber E., Shekhtman R., 2013. Three dimensional inversion of multisource time domain electromagnetic data, Geophysics , 78, E47– E57. https://doi.org/10.1190/geo2012-0131.1 Google Scholar CrossRef Search ADS   Persova M.G., Soloveichik Y.G., Trigubovigh G.M., 2011. Computer modeling of geoelectromagnetic fields in three-dimensional media by the finite element method, Izv. Phys. Solid Earth , 47, 79– 89. https://doi.org/10.1134/S1069351311010095 Google Scholar CrossRef Search ADS   Persova M.G., Soloveichik Y.G., Trigubovich G.M., Tokareva M.G., 2013. Methods and algorithms for reconstructing three-dimensional distributions of electric conductivity and polarization in the medium by finite-element 3D modeling using the data of electromagnetic sounding, Izv. Phys. Solid Earth , 49, 329– 343. https://doi.org/10.1134/S1069351313030117 Google Scholar CrossRef Search ADS   Persova M.G., Soloveichik Y.G., Domnikov P.A., Vagin D.V., 2014a. Finite element modeling of three-dimensional geoelectromagnetic fields excited by arbitrarily oriented induction source, in 12th International Conference on Actual Problems of Electronic Instrument Engineering, APEIE 2014 , pp. 598– 602, Institute of Electrical and Electronics Engineers Inc. Persova M.G., Soloveichik Y.G., Trigubovich G.M., Vagin D.V., Domnikov P.A., 2014b. Transient electromagnetic modelling of an isolated wire loop over a conductive medium, Geophys. Prospect. , 62, 1193– 1201. https://doi.org/10.1111/1365-2478.12122 Google Scholar CrossRef Search ADS   Persova M.G., Soloveichik Y.G., Domnikov P.A., Vagin D.V., Koshkina Y.I., 2015a. Electromagnetic field analysis in the marine CSEM detection of homogeneous and inhomogeneous hydrocarbon 3D reservoirs, J. Appl. Geophys. , 119, 147– 155. https://doi.org/10.1016/j.jappgeo.2015.05.019 Google Scholar CrossRef Search ADS   Persova M.G., Soloveichik Y.G., Simon E.I., Koshkina Y.I., Epanchintseva T.B., 2015b. Methods and software to perform 3D-inversion of the airborne electrical prospecting data in time domain, in Geophysics 2015—11th EAGE International Scientific and Practical Conference and Exhibition on Engineering and Mining Geophysics , doi:10.3997/2214-4609.201412271. Persova M.G., Soloveichik Y.G., Koshkina Y.I., Vagin D.V., Trubacheva O.S., 2016. Geometrical Nonlinear 3D Inversion of Airborne Time Domain EM Data, in Near Surface Geoscience 2016—First Conference on Geophysics for Mineral Exploration and Mining , Barcelona, Spain. Puzyrev V., Koric S., Wilkin S., 2016. Evaluation of parallel direct sparse linear solvers in electromagnetic geophysical problems, Comput. Geosci. , 89, 79– 87. https://doi.org/10.1016/j.cageo.2016.01.009 Google Scholar CrossRef Search ADS   Schenk O., G̈artner K., 2004. Solving unsymmetric sparse systems of linear equations with PARDISO, Future Gener. Comput. Syst. , 20, 475– 487. https://doi.org/10.1016/j.future.2003.07.011 Google Scholar CrossRef Search ADS   Schwarzbach C., Börner R.-U., Spitzer K., 2011. Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics - a marine CSEM example, Geophys. J. Int. , 187, 63– 74. https://doi.org/10.1111/j.1365-246X.2011.05127.x Google Scholar CrossRef Search ADS   Schwarzbach C., Haber E., 2013. Finite element based inversion for time-harmonic electromagnetic problems, Geophys. J. Int. , 193, 615– 634. https://doi.org/10.1093/gji/ggt006 Google Scholar CrossRef Search ADS   Shantsev D.V., Maaø F.A., 2015. Rigorous interpolation near tilted interfaces in 3-D finite-difference EM modelling, Geophys. J. Int. , 200, 743– 755. https://doi.org/10.1093/gji/ggu429 Google Scholar CrossRef Search ADS   Soloveichik Y.G., Royak M.E., Moiseev V.S., Vasilev A.V., 1997. Finite element modeling of 3D electric fields in electrical prospecting, Fiz. Zemli , 9, 67– 71. Soloveichik Y.G., Royak M.E., Moiseev V.S., Trigubovich G.M., 1998. Three-dimensional modeling of nonstationary electromagnetic fields using the finite element method, Fiz. Zemli , 10, 78– 83. Soloveichik Y.G., Persova M.G., Domnikov P.A., Vagin D.V., 2014. Method for calculating the three-dimensional time-harmonic electromagnetic fields in marine electrical prospecting, in 12th International Conference on Actual Problems of Electronic Instrument Engineering, APEIE 2014 , pp. 594– 597, Institute of Electrical and Electronics Engineers Inc. Ullmann A., Scheunert M., Afanasjew M., Börner R.U., Siemon B., Spitzer K., 2016. A cut-&-paste strategy for the 3-D inversion of helicopter-borne electromagnetic data—II. Combining regional 1-D and local 3-D inversion, J. Appl. Geophys. , 130, 131– 144. https://doi.org/10.1016/j.jappgeo.2016.04.008 Google Scholar CrossRef Search ADS   Um E.S., Commer M., Newman G.A., Hoversten G.M., 2015. Finite element modelling of transient electromagnetic fields near steel-cased wells, Geophys. J. Int. , 202, 901– 913. https://doi.org/10.1093/gji/ggv193 Google Scholar CrossRef Search ADS   Usui Y., 2015. 3-D inversion of magnetotelluric data using unstructured tetrahedral elements: applicability to data affected by topography, Geophys. J. Int. , 202, 828– 849. https://doi.org/10.1093/gji/ggv186 Google Scholar CrossRef Search ADS   Wang W., Wu X., Spitzer K., 2013. Three-dimensional DC anisotropic resistivity modelling using finite elements on unstructured grids, Geophys. J. Int. , 193, 734– 746. https://doi.org/10.1093/gji/ggs124 Google Scholar CrossRef Search ADS   Yang D., Oldenburg D.W., Haber E., 2014. 3-D inversion of airborne electromagnetic data parallelized and accelerated by local mesh and adaptive soundings, Geophys. J. Int. , 196, 1492– 1507. https://doi.org/10.1093/gji/ggt465 Google Scholar CrossRef Search ADS   Yang K., Yılmaz A.E., Torres-Verdín C., 2017. A goal-oriented framework for rapid integral-equation-based simulation of borehole resistivity measurements of 3D hydraulic fractures, Geophysics , 82, D123– D133. https://doi.org/10.1190/geo2016-0179.1 Google Scholar CrossRef Search ADS   Yin C., Zhang B., Liu Y., Cai J., 2016. A goal-oriented adaptive finite-element method for 3D scattered airborne electromagnetic method modeling, Geophysics , 81, E337– E346. https://doi.org/10.1190/geo2015-0580.1 Google Scholar CrossRef Search ADS   Zaslavsky M., Druskin V., Davydycheva S., Knizhnerman L., Abubakar A., Habashy T., 2011. Hybrid finite-difference integral equation solver for 3D frequency domain anisotropic electromagnetic problems, Geophysics , 76, F123– F137. https://doi.org/10.1190/1.3552595 Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

Finite-element solution to multidimensional multisource electromagnetic problems in the frequency domain using non-conforming meshes

Loading next page...
 
/lp/ou_press/finite-element-solution-to-multidimensional-multisource-vNc37b80J0
Publisher
The Royal Astronomical Society
Copyright
© The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx530
Publisher site
See Article on Publisher Site

Abstract

Summary We propose an approach to solving multisource induction logging problems in multidimensional media. According to the type of induction logging tools, the measurements are performed in the frequency range of 10 kHz to 14 MHz, transmitter–receiver offsets vary in the range of 0.5–8 m or more, and the trajectory length is up to 1 km. For calculating the total field, the primary-secondary field approach is used. The secondary field is calculated with the use of the finite-element method (FEM), irregular non-conforming meshes with local refinements and a direct solver. The approach to constructing basis functions with the continuous tangential components (from Hcurl(Ω)) on the non-conforming meshes from the standard shape vector functions is developed. On the basis of this method, the algorithm of generating global matrices and a vector of the finite-element equation system is proposed. We also propose the method of grouping the logging tool positions, which makes it possible to significantly increase the computational effectiveness. This is achieved due to the compromise between the possibility of using the 1-D background medium, which is very similar to the investigated multidimensional medium for a small group, and the decrease in the number of the finite-element matrix factorizations with the increasing number of tool positions in one group. For calculating the primary field, we propose the method based on the use of FEM. This method is highly effective when the 1-D field is required to be calculated at a great number of points. The use of this method significantly increases the effectiveness of the primary-secondary field approach. The proposed approach makes it possible to perform modelling both in the 2.5-D case (i.e. without taking into account a borehole and/or invasion zone effect) and the 3-D case (i.e. for models with a borehole and invasion zone). The accuracy of numerical results obtained with the use of the proposed approach is compared with the one obtained by other codes for 1-D and 3-D anisotropic models. The results of this comparison lend support to the validity of our code. We also present the numerical results proving greater effectiveness of the finite-element approach proposed for calculating the 1-D field in comparison with the known codes implementing the semi-analytical methods for the case in which the field is calculated at a large number of points. Additionally, we present the numerical results which confirm the accuracy advantages of the automatic choice of a background medium for calculating the 1-D field as well as the results of 2.5-D modelling for a geoelectrical model with anisotropic layers, a fault and long tool-movement trajectory with the varying dip angle. Electrical properties, Downhole methods, Numerical approximations and analysis, Numerical modelling 1 INTRODUCTION For processing the electromagnetic soundings data, the multidimensional modelling is increasingly used (e.g. Oldenburg et al.2013; Kamm & Pedersen 2014; Kordy et al.2016b). The inversions based on the 2-D or 3-D modelling make it possible to recover the structure of the investigated medium with greater adequacy (e.g. Persova et al.2013; Schwarzbach & Haber 2013; Haber & Schwarzbach 2014; Usui 2015; Persova et al.2016, etc.). However, these inversions require high computational costs. Therefore in practice, until recently, quite simple and low-cost approaches, though not always adequate, have often been used. They are based on calculating the electromagnetic fields in the 1-D medium or other approaches, which make use of simplified mathematical models of the electromagnetic field (among others Macnae 2015; Hunkeler et al.2016). In some cases, these approaches allow us to obtain acceptable results (e.g. Hennessy & Macnae 2010; Hunkeler et al.2016). Nevertheless, the researchers often note that the simplified and 1-D approaches not rarely give the distorted results, especially near local inhomogeneities (e.g. Ley-Cooper et al.2015; Persova et al.2015b, Ullmann et al.2016). In induction logging, this problem is important, e.g. when investigating the horizontal wells, and if the geometry and electrophysical properties of the borehole environment are locally and drastically changed or it is necessary to correctly define the intersection of the fault by the borehole. In these cases, for adequately recovering the medium, it is necessary to use the multidimensional approaches based on accurate (without significant simplifications) mathematical models of the electromagnetic field. Therefore, the increasingly important problem is the development of such methods of electromagnetic fields calculation that, on the one hand, would make it possible to model the field in a complex medium with the required adequacy (with the use of mathematical models without simplifications) and, on the other hand, would be computationally low-cost and convenient to use. For multidimensional electromagnetic fields modelling in complex 3-D media, among the mesh-based numerical methods the following are more frequently used: the finite-difference (FD) method (Grayver et al.2013; Jaysaval et al.2015; Shantsev & Maaø 2015, etc.), the finite volume method (Haber & Ruthotto 2014; Haber & Schwarzbach 2014; Jahandari & Farquharson 2014; Jahandari & Farquharson 2015, etc.) and the finite-element (FE) method (Farquharson & Miensopust 2011; Ansari & Farquharson 2014; Börner et al.2015; Persova et al.2015a; Kordy et al.2016a, etc.). Besides, for modelling the multidimensional electromagnetic fields including the solution of 3-D induction logging problems, the integral equation method is applied quite successfully (Fang et al.2003; Nie et al.2013; Yang et al.2017, etc.). Its computational effectiveness is especially high if the volume of 3-D inhomogeneities is small. When solving the multidimensional electromagnetic problems with the use of mesh-based numerical methods, a significant lowering of computational costs can be achieved by means of the following approaches. First, this is the primary-secondary field approach which uses the separate calculation of the horizontally layered medium field (that is called the primary field). Among the first papers, which considered this approach for FD and FE methods, one should mention the papers (Newman & Alumbaugh 1995; Soloveichik et al.1997; Soloveichik et al.1998; Badea et al.2001). Nowadays, this approach continues to be actively applied (e.g. Mukherjee & Everett 2011; Zaslavsky et al.2011; Da Silva et al.2012; Grayver et al.2013). The more similar the investigated medium to the horizontally layered one, which can be taken different for different logging tool positions, the higher the effectiveness of the primary-secondary field approach. In this case, for achieving the required accuracy of the summary (total) field calculation, when the anomalous (secondary) field is calculated, the number of nodes in the 3-D mesh can be decreased by several times, and the computational time of one 3-D problem solution can be decreased by an order or more. Second, a significant lowering of the computational costs can be achieved if we use the unstructured meshes (e.g. Lelièvre & Farquharson 2013; Wang et al.2013; Um et al.2015) including the non-conforming ones with local refinements (e.g. Grayver & Bürg 2014; Persova et al.2014b; Grayver 2015). In this case, in comparison with the regular meshes, the number of elements in the mesh can be decreased by several times (without the loss of the accuracy of the 3-D problem solution), and this allows the additional reduction of the computational costs by an order. Third, if the number of the logging tool positions is big, a significant lowering of the computational costs can be achieved if we use the direct solvers for the finite-difference or FE equations system (Da Silva et al.2012; Oldenburg et al.2013; Chung et al.2014; Haber & Schwarzbach 2014; Puzyrev et al.2016) (e.g. PARDISO solver (Schenk & G̈artner 2004)). In this case, the matrix factorization, which is the most computational cost procedure in calculating a 3-D field, is performed only once for several logging tool positions simultaneously, and the computational time of 3-D problems for all logging tool positions joined in one relatively small group is close to the computational time of the 3-D problem for one tool position. However, it is necessary to take into account that the computational time for one tool position when the direct solver is used can be much higher in comparison with the time when the iteration solver is used. In this paper, we consider the approach to calculating the 3-D field in the induction logging problems, in which all of the three above methods are used for lowering the computational costs. Additionally, we propose: a new method of FE calculation of the primary field in a horizontally layered medium, which is effective when it is necessary to calculate the primary field at many points simultaneously (that is necessary for calculating the secondary field); the method of 3-D electromagnetic field FE calculation on irregular non-conforming meshes with the construction of conforming basis functions from the non-conforming ones and corresponding procedure of assembling matrices of FEs; the method of grouping the logging tool positions (in order to use the direct solver when calculating the 3-D field for all positions in one group) with the automatic choice of the 1-D background medium for calculating the primary fields and automatic generation of the irregular 3-D mesh with hanging nodes for one group of logging tool positions. 2 MATHEMATICAL MODELS 2.1 The model describing the electromagnetic field in a 3-D anisotropic medium We represent a total electric field in the form $${{{\bf \vec{E}}}^t} = {{{\bf \vec{E}}}^p} + {{{\bf \vec{E}}}^s}$$ where $${{{\bf \vec{E}}}^p}$$ is the primary electric field induced by a source in a horizontally layered (1-D background) medium, and $${{{\bf \vec{E}}}^s}$$ is the 3-D secondary electric field determined by the impact of 3-D inhomogeneities such as layers pinchout, faults, non-parallelism of layers, etc. We can seek $${{{\bf \vec{E}}}^s}$$ by solving the following vector differential equation:   \begin{equation} \frac{1}{{{\mu _0}}}\nabla \times \nabla \times {{{\bf \vec{E}}}^s} + i\omega \sigma {{{\bf \vec{E}}}^s} - {\omega ^2}\varepsilon {{{\bf \vec{E}}}^s} = - i\omega \left( {\sigma - {\sigma ^p}} \right){{{\bf \vec{E}}}^p} + {\omega ^2}\varepsilon {{{\bf \vec{E}}}^p}, \end{equation} (1)where μ0 denotes the vacuum magnetic permeability, σ denotes the 3 × 3 conductivity tensor in a 3-D medium, σp denotes the 3 × 3 conductivity tensor in a horizontally layered (1-D background) medium (for which $${{{\bf \vec{E}}}^p}$$ is calculated), ε(x, y, z) denotes the dielectric permeability of the 3-D medium, ω = 2πν denotes the angular frequency (ν is the electromagnetic field frequency), and $$i = \sqrt { - 1} $$. The components σij(x, y, z) of the tensor σ and components $$\sigma _{ii}^p( z )$$ of the tensor σp are the piecewise constant functions of coordinates. The components σii are not equal to $$\sigma _{ii}^p$$ only in the subdomains where 3-D inhomogeneities are located. When i ≠ j, σij are not equal to zero only in the subdomains where 3-D anisotropic inhomogeneities are located, and the anisotropic axes of these inhomogeneities do not coincide with the coordinate axes. Dielectric permeability ε can also be defined as a tensor. If thin layers in the anisotropic medium are not parallel to Cartesian planes (i.e. anisotropic axes do not coincide with coordinate axes), the conductivity tensor can be obtained by the relation   \begin{equation} \sigma = {{{\bf M}}^{{\bf T}}}\tilde{\sigma }{{\bf M}}, \end{equation} (2)where M is the rotation matrix, $$\tilde{\sigma }$$ is the diagonal tensor with the components corresponding to the lateral and vertical conductivities (in the local coordinate system related with anisotropic axes), the superscript ‘T’ denotes transposition. In a similar way, the tensor of dielectric permeability can be defined in the anisotropic medium. Eq. (1) defines the secondary field for the case when the primary field is calculated without the account of displacement currents. In general, in the induction logging problems, when electromagnetic field is calculated for the frequency ν < 200 kHz, the displacement currents effects are minor, and the terms with ε can be excluded from eq. (1) (i.e. ε is assumed to be zero). If the frequency ν > 200 kHz, the displacement currents effect on the calculated electromagnetic field can be significant. In this case, the displacement currents are to be taken into account in calculating the secondary field. In induction logging measurements, the transmitters (induction coils) with different orientation can be used. Most often the coils are coaxial with the borehole. But sometimes the tools with several (two or three) orthogonal coils (the axis of one of these coils coincides with the tool axis) are considered. For computing $${{{\bf \vec{E}}}^p}$$ induced by an arbitrary oriented coil in the horizontally layered background medium, it is required to compute two fields for two coil orientations: in the first case the coil axis is perpendicular to the layers, in the second one the coil axis is parallel to the layers. In general, for calculating $${{{\bf \vec{E}}}^p}$$, analytical or semi-analytical methods can be used. The first and the simplest is the approach where the field of homogeneous space is taken as $${{{\bf \vec{E}}}^p}$$. This field can be calculated by the analytical method. The disadvantage of this approach is the following. Even if the conductivity of the homogeneous space is taken equal to the conductivity of the layer in which the transmitter is located, it is impossible to reach a comparable effectiveness that can be achieved if we take the horizontally layered medium field as $${{{\bf \vec{E}}}^p}$$. This is because, in the first case, the layers (as 3-D inhomogeneities) are the sources of the secondary field. Therefore, for reaching the required accuracy, it is necessary to use sufficiently fine meshes (that will lead to much higher computational costs). Until recently, semi-analytical methods have been developed quite intensively (Key 2009). They are very effective if there is a need to compute the field at a relatively moderate number of points (corresponding to receiver locations). But if it is required to compute the field values at a great number of points, the computational efficiency of these methods is not so high. For implementing the method of calculating the 3-D electromagnetic field when the horizontally layered medium is taken as a background medium, it is required to compute the primary field $${{{\bf \vec{E}}}^p}$$ at a great number of points (cells) in subdomains (of the computational domain) where 3-D inhomogeneities are located. Even if the optimized irregular non-conforming meshes are used, the primary field needs to be calculated at tens of thousands of points. Consequently, in this case the approaches which allow us to compute the horizontally layered medium field in all the investigated volume simultaneously are more advantageous (Persova et al.2011; Mogilatov et al.2016). Therefore, for calculating the primary field $${{{\bf \vec{E}}}^p}$$, we propose another approach based on the use of the FE method. 2.2 Calculating $${{{\bf \vec{E}}}^p}$$ in a horizontally layered background medium We represent $${{{\bf \vec{E}}}^p}$$ as a sum of solutions of several axisymmetrical problems. When describing mathematical models in the cylindrical coordinate system, we use the diagonal 2 × 2 tensor σp, rz(z) with the components $$\sigma _r^{p,rz} = \sigma _{11}^p = \sigma _{22}^p$$ and $$\sigma _z^{p,rz} = \sigma _{33}^p$$ (where $$\sigma _{ii}^p$$ are the diagonal components of the tensor σp in eq. 1). At first, we consider the mathematical model for computing the electrical field $${{{\bf \vec{E}}}^{p,1}}$$ which is induced by a ‘horizontal’ coil. This coil lies in plane which is parallel to the layers of a horizontally layered background medium (and, obviously, the coil axis is perpendicular to the layers borders). $${{{\bf \vec{E}}}^{p,1}}$$ induced by this coil is axisymmetrical and has only one non-zero component $$E_\varphi ^{}( {r,z} )$$ in cylindrical coordinates {r, φ, z}. This component ($$E_\varphi ^{}( {r,z} )$$) can be found by solving the equation   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta E_\varphi ^{} + \frac{1}{{{\mu _0}{r^2}}}E_\varphi ^{} + i\omega \sigma _r^{p,rz}E_\varphi ^{} = - i\omega I{\delta ^{{J_\varphi }}}, \end{equation} (3)where I is the value of the coil current, and $${\delta ^{{J_\varphi }}}$$ is the delta-function concentrated on the circle r = rδ, z = zδ corresponding to the coil contour ($$I{\delta ^{{J_\varphi }}}$$ is the source current density Jφ). The components of the vectors $${{{\bf \vec{E}}}^{p,1}}$$ and $${{{\bf \vec{B}}}^{p,1}}$$ can be obtained using $$E_\varphi ^{}$$ by the following relations:   \begin{equation} E_x^{p,1}\left( {x,y,z} \right) = - E_\varphi ^{}\left( {r,z} \right)\frac{{y - {y^c}}}{r},\quad E_y^{p,1}\left( {x,y,z} \right) = E_\varphi ^{}\left( {r,z} \right)\frac{{x - {x^c}}}{r},\quad E_z^p = 0, \end{equation} (4)  \begin{equation*} B_x^{p,1}\left( {x,y,z} \right) = \frac{1}{{i\omega }}\frac{{\partial {E_\varphi }\left( {r,z} \right)}}{{\partial z}}\frac{{x - {x^c}}}{r},\quad B_y^{p,1}\left( {x,y,z} \right) = \frac{1}{{i\omega }}\frac{{\partial {E_\varphi }\left( {r,z} \right)}}{{\partial z}}\frac{{y - {y^c}}}{r}, \end{equation*}   \begin{equation} B_z^{p,1}\left( {x,y,z} \right) = - \frac{1}{{i\omega }}\left( {\frac{{{E_\varphi }\left( {r,z} \right)}}{r} + \frac{{\partial {E_\varphi }\left( {r,z} \right)}}{{\partial r}}} \right), \end{equation} (5)where xc and yc are the coordinates of the coil centre, and $$r = \sqrt {{{( {x - {x^c}} )}^2} + {{( {y - {y^c}} )}^2}} $$. The electrical field $${{{\bf \vec{E}}}^{p,2}}$$ induced by a ‘vertical’ coil (its axis is parallel to the borders of the layers of the horizontally layered medium) is 3-D even for a 1-D (horizontally layered) medium, that is, the components of $${{{\bf \vec{E}}}^{p,2}}$$ depend on all three space coordinates. In order to compute this field, we apply the approach similar to the one considered by Persova et al. (2011) for modelling the field induced by the horizontal electrical dipole. It is based on the representation of the field, which is in fact a 3-D field in the 1-D medium, as a sum of the fields each of which can be obtained by the solution of the problem of a less dimension. However, when we use the approach considered by Persova et al. (2011) there appears a problem (this is particularly true for harmonic modes) which is related to the fact that the total field is formed by the components with the opposite signs the moduli of which are close to each other. This makes it necessary to increase significantly the accuracy of the components calculation that, in its turn, leads to the significant increase in computational costs. Therefore, in this paper we propose the modification of the approach described by Persova et al. (2011) that makes it possible to obtain $${{{\bf \vec{E}}}^{p,2}}$$ (the field induced by the ‘vertical’ coil) as a sum of several axisymmetrical fields without mutually subtracting components (Soloveichik et al.2014; Persova et al.2014a). Each of these fields (forming $${{{\bf \vec{E}}}^{p,2}}$$) is induced by a special source and can be obtained by solving the corresponding axisymmetrical problem. At first we note that, in the problem considered, the sizes of the transmitter coil are small relative to the distance to the receivers. Therefore, we can take any geometrical form of coil. We represent the ‘vertical’ coil as a rectangular (square) frame. In this case the total source (frame) can be represented as a sum of sources shown in Fig. 1. This representation allows us to obtain the electromagnetic field in a horizontally layered background medium by means of solutions of several axisymmetrical problems. Figure 1. View largeDownload slide The view of (a) non-axisymmetric source ‘I1’ and (b) two axisymmetrical sources ‘I2’; (c) the current of the upper part of the source ‘I1’ in plane XY; (d) the currents of one of NGR sources in plane XY; (e) currents configuration which is obtained when the negative source ‘I2’ is added to the source ‘I1’; (f) currents configuration which is obtained when the positive source ‘I2’ is added to the source shown in panel (e) (as a result, the ‘vertical’ rectangular current frame is obtained); (g) the approximation of the horizontal ‘current lines’ by the point sources located in the centres of subintervals of these ‘current lines’ (the length of subintervals is $$h_L^{}$$). Dashed light-green arrows show the direction of the surface radial currents flowing in plane XY. Solid black lines show the currents flowing in the ‘vertical’ current frame. Red points show the positions of the point sources approximating the horizontal ‘current lines’. Figure 1. View largeDownload slide The view of (a) non-axisymmetric source ‘I1’ and (b) two axisymmetrical sources ‘I2’; (c) the current of the upper part of the source ‘I1’ in plane XY; (d) the currents of one of NGR sources in plane XY; (e) currents configuration which is obtained when the negative source ‘I2’ is added to the source ‘I1’; (f) currents configuration which is obtained when the positive source ‘I2’ is added to the source shown in panel (e) (as a result, the ‘vertical’ rectangular current frame is obtained); (g) the approximation of the horizontal ‘current lines’ by the point sources located in the centres of subintervals of these ‘current lines’ (the length of subintervals is $$h_L^{}$$). Dashed light-green arrows show the direction of the surface radial currents flowing in plane XY. Solid black lines show the currents flowing in the ‘vertical’ current frame. Red points show the positions of the point sources approximating the horizontal ‘current lines’. The rectangular ‘vertical’ current frame is represented as a sum of source ‘I1’ shown in Fig. 1(a) (we consider the top and bottom currents as the parts of one source ‘I1’) and two sources ‘I2’ shown in Fig. 1(b). The source ‘I1’ consists of four surface radial currents, which flow in plane XY (in Fig. 1, their direction is shown by dashed light-green arrows), and two horizontal ‘current lines’ (in Fig. 1, the currents in them are shown by solid black lines). We consider the left and right currents as different sources of the type ‘I2’, the former is negative and the latter is positive with regard to the direction of the axes shown in Fig. 1 (in this case we direct the z-axis down because in all examples presented in this paper, the direction of the coordinate axes is the same). Each of the sources ‘I2’ consists of the vertical ‘current line’ (the current in which shown by solid black line) and two surface radial currents (their direction is shown by dashed light-green arrows). Joining the top and bottom currents in one source ‘I1’ and separating the left and right currents into two sources of the type ‘I2’ are convenient for further reduction of the fields induced by these sources to the corresponding axisymmetrical fields. The Z-coordinate of the upper part of source ‘I1’ coincides with the Z-coordinate of the upper frame border, and the Z-coordinate of the lower part of the source ‘I1’ (with the opposite sign of the current) coincides with the Z-coordinate of the lower frame border. The axis of one of the sources ‘I2’ coincides with the left frame border, and the axis of another source ‘I2’ (with an opposite sign of the current) coincides with the right frame border. Convergent and divergent surface radial currents are the parts of these sources. They provide zero current divergence in sources ‘I1’ and ‘I2’. Further, we will refer to as ‘non-grounded radial sources’ (NGR sources) for each of these surface radial currents converging into the point Pi or diverging from it. Fig. 1(c) shows the current of the upper part of the source ‘I1’ in plane XY, and Fig. 1(d) shows the currents of one of NGR sources in plane XY. For each of NGR sources, the same current passes through any circle C the centre of which coincides with the NGR source centre (see Fig. 1d, the integrals along circles С1 and С2 of surface radial currents of the NGR source are equal). Figs 1(e) and (f) demonstrate how the current in the ‘vertical’ frame is composed from the source ‘I1’ and two sources ‘I2’. In Fig. 1(e) the negative source ‘I2’ with the axis, which coincides with the left side of the ‘vertical’ frame, is added to the source ‘I1’. The radial currents of the source ‘I2’ compensate the radial currents of the source ‘I1’ in the left top and left bottom angles (these currents are equal and are in the opposite directions). When the positive source ‘I2’ (with the axis coinciding with the right side of the ‘vertical’ frame) is added to the source shown in Fig. 1(e), the radial currents in the right top and right bottom angles are compensated. As a result, we obtain the source with the current in the ‘vertical’ frame shown in Fig. 1(f). We would like to reiterate that each of the sources ‘I2’ generates the axisymmetrical field (in the coordinate system in which the z-axis coincides with the vertical current of the source). The mathematical model for calculating the axisymmetrical field generated by the source ‘I2’ is presented below (see eq. 21–23). We first consider how the field induced by the non-axisymmetric source ‘I1’ can be calculated by solving several axisymmetrical problems. We introduce the vector potential $${{\bf \vec{A}}}$$ that is related with the magnetic induction and electric field by the relationships $${{\bf \vec{B}}} = {\rm{rot}}{{\bf \vec{A}}}$$ and $${{\bf \vec{E}}} = - \frac{{\partial {{\bf \vec{A}}}}}{{\partial t}} = - i\omega {{\bf \vec{A}}}$$. Then, the vector potential $${{{\bf \vec{A}}}^{{\rm{I1}}}}$$ defining the electromagnetic field induced by source ‘I1’ in a horizontally layered medium (with conductivity σp) can be found by solving the equation   \begin{equation} \frac{1}{{{\mu _0}}}\nabla \times \nabla \times {{{\bf \vec{A}}}^{{\rm{I1}}}} + i\omega {\sigma ^p}{{{\bf \vec{A}}}^{{\rm{I1}}}} = {{{\bf \vec{J}}}^{{\rm{I1}}}}, \end{equation} (6)where $${{{\bf \vec{J}}}^{{\rm{I1}}}}$$ is the density of the currents defining source ‘I1’. We define this source in such a way that, first, its currents have only horizontal non-zero components, and, second, for any bounded subdomain Ωf, the value of the current entering Ωf is equal to the value of the outward current. Therefore, for the current density $${{{\bf \vec{J}}}^{{\rm{I1}}}}$$, the ratio   \begin{equation} \!\nabla \cdot {{{\bf \vec{J}}}^{{\rm{I1}}}} = 0 \end{equation} (7)is true and $$J_z^{{\rm{I1}}} \equiv 0$$. It would now be easy to show that eq. (6) is satisfied for the vector potential $${{{\bf \vec{A}}}^{{\rm{I1}}}}$$ whose z-component is equal to zero and the components $$A_x^{{\rm{I1}}}$$ and $$A_y^{{\rm{I1}}}$$ satisfy the equations   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta A_x^{{\rm{I1}}} + i\omega \sigma _{11}^pA_x^{{\rm{I1}}} = J_x^{{\rm{I1}}}, \end{equation} (8)  \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta A_y^{{\rm{I1}}} + i\omega \sigma _{22}^pA_y^{{\rm{I1}}} = J_y^{{\rm{I1}}}, \end{equation} (9)where Δ denotes the Laplacian, and $$\sigma _{11}^p = \sigma _{22}^p = \sigma _r^{p,rz}(z)$$. We will prove it. For this purpose, we differentiate eq. (8) with respect to x, and eq. (9) with respect to y (taking into account that $$\sigma _{11}^p$$ and $$\sigma _{22}^p$$ do not depend on x and y and equal to $$\sigma _r^{p,rz}$$). Changing the order of differentiation in the first summands of the left-hand sides of these equations and summing them, we obtain   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta \left( {\frac{{\partial\! A_x^{{\rm{I1}}}}}{{\partial x}} + \frac{{\partial\! A_y^{{\rm{I1}}}}}{{\partial y}}} \right) + i\omega \sigma _r^{p,rz}\left( {\frac{{\partial\! A_x^{{\rm{I1}}}}}{{\partial x}} + \frac{{\partial\! A_y^{{\rm{I1}}}}}{{\partial y}}} \right) = \left( {\frac{{\partial J_x^{{\rm{I1}}}}}{{\partial x}} + \frac{{\partial J_y^{{\rm{I1}}}}}{{\partial y}}} \right). \end{equation} (10)Taking into account that $$A_z^{{\rm{I1}}} \equiv 0$$ and $$J_z^{{\rm{I1}}} \equiv 0$$ (and, consequently, $$\frac{{\partial\! A_z^{{\rm{I1}}}}}{{\partial z}} = 0$$ and $$\frac{{\partial J_z^{{\rm{I1}}}}}{{\partial z}} = 0$$), we obtain   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta ( {\nabla \cdot {{{{\bf \vec{A}}}}^{{\rm{I1}}}}}) + i\omega \sigma _r^{p,rz}\left( {\nabla \cdot {{{{\bf \vec{A}}}}^{{\rm{I1}}}}} \right) = \nabla \cdot {{{\bf \vec{J}}}^{{\rm{I1}}}}. \end{equation} (11)As far as eq. (7) is true for $${{{\bf \vec{J}}}^{{\rm{I1}}}}$$, the function $$u = \nabla \cdot {{{\bf \vec{A}}}^{{\rm{I1}}}}$$ satisfies the homogeneous equation   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta u + i\omega \sigma _r^{p,rz}u = 0 \end{equation} (12)for which the function u = 0 is the solution. Thus, the divergence of vector potential $${{{\bf \vec{A}}}^{{\rm{I1}}}}$$ (whose x- and y-components satisfy eqs (8) and (9), and $$A_z^{{\rm{I1}}} \equiv 0$$) is equal to zero. Taking into account $$A_z^{{\rm{I1}}} = 0$$, the equation system (8), (9) can be written in the vector form   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta {{{\bf \vec{A}}}^{{\rm{I1}}}} + i\omega {\sigma ^p}{{{\bf \vec{A}}}^{{\rm{I1}}}} = {{{\bf \vec{J}}}^{{\rm{I1}}}}, \end{equation} (13)where Δ denotes the vector Laplacian. As far as $$\nabla \cdot {{{\bf \vec{A}}}^{{\rm{I1}}}} = 0$$ (we have just shown it), and taking into account the identical relation $$\nabla \times \nabla \times {{{\bf \vec{A}}}^{{\rm{I1}}}} \equiv - \Delta {{{\bf \vec{A}}}^{{\rm{I1}}}} + \nabla \nabla \cdot {{{\bf \vec{A}}}^{{\rm{I1}}}}$$, the vector potential $${{{\bf \vec{A}}}^{{\rm{I1}}}}$$, which is obtained by solving eqs (8) and (9) (or by solving vector eq. 13), satisfies eq. (6). Now we show how the solution of eq. (13) could be obtained by solving several axisymmetrical problems. For this purpose, we divide the source ‘I1’ into the following parts: four NGR sources and two ‘current lines’. By ‘current line’ we mean the source which is concentrated on one of the horizontal borders P1P2 or P3P4 of the ‘vertical’ current frame (see Fig. 1). The density of such sources can be defined via delta-function δL, which is concentrated on the line and defined in the following way: the integral of δL taken over an arbitrary volume Ωf is equal to the length of the L line segment, which lies inside the volume Ωf. Thus, we represent right-hand side of vector eq. (13) as a sum of three summands:   \begin{equation} {{{\bf \vec{J}}}^{{\rm{I1}}}} = {{{\bf \vec{J}}}^{{\rm{I1,NGR1}}}} + {{{\bf \vec{J}}}^{{\rm{I1,NGR2}}}} + {{{\bf \vec{J}}}^{{\rm{I1,hor}}}}, \end{equation} (14)where $${{{\bf \vec{J}}}^{{\rm{I1,NGR1}}}}$$ is the density of surface currents in two (axial) NGR sources, while the currents of one of these sources converge into the point P1, and currents of another source diverge from the point P3 (see Fig. 1a); $${{{\bf \vec{J}}}^{{\rm{I1,NGR2}}}}$$ is the density of surface currents in other two (axial) NGR sources (currents of one of these sources diverge from the point P2, and currents of another source converge to the point P4); $${{{\bf \vec{J}}}^{{\rm{I1,hor}}}} = {{{\bf \vec{J}}}^{{\rm{I1,}}{{\rm{P}}_1}{{\rm{P}}_2}}} - {{{\bf \vec{J}}}^{{\rm{I1,}}{{\rm{P}}_3}{{\rm{P}}_4}}}$$ is the density of currents in the horizontal ‘current lines’ P1P2 and P3P4. Then, the solution of eq. (13) can be found as a sum of solutions of three equations (they are in fact independent)   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta {{{\bf \vec{A}}}^{{\rm{I1,NGR1}}}} + i\omega {\sigma ^p}{{{\bf \vec{A}}}^{{\rm{I1,NGR1}}}} = {{{\bf \vec{J}}}^{{\rm{I1,NGR1}}}}, \end{equation} (15)  \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta {{{\bf \vec{A}}}^{{\rm{I1,NGR2}}}} + i\omega {\sigma ^p}{{{\bf \vec{A}}}^{{\rm{I1,NGR2}}}} = {{{\bf \vec{J}}}^{{\rm{I1,NGR2}}}}, \end{equation} (16)  \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta {{{\bf \vec{A}}}^{{\rm{I1,hor}}}} + i\omega {\sigma ^p}{{{\bf \vec{A}}}^{{\rm{I1,hor}}}} = {{{\bf \vec{J}}}^{{\rm{I1,}}{{\rm{P}}_1}{{\rm{P}}_2}}} - {{{\bf \vec{J}}}^{{\rm{I1,}}{{\rm{P}}_3}{{\rm{P}}_4}}}. \end{equation} (17)Right-hand sides of eqs (15) and (16) are the functions with axial symmetry (the z-axis passes through the points P1,P3 or P2,P4, respectively). In the corresponding cylindrical coordinate system, the functions $${{{\bf \vec{J}}}^{{\rm{I1,NGR1}}}}$$ and $${{{\bf \vec{J}}}^{{\rm{I1,NGR2}}}}$$ of the right-hand sides of these equations have the only non-zero radial component: $$J_r^{{\rm{I1,NGR1}}} = \frac{I}{{2\pi r}}( {{\delta ^{{{\hat{\Gamma }}_2}}} - {\delta ^{{{\hat{\Gamma }}_1}}}} )$$ and $$J_r^{{\rm{I1,NGR2}}} = \frac{I}{{2\pi r}}( {{\delta ^{{{\hat{\Gamma }}_1}}} - {\delta ^{{{\hat{\Gamma }}_2}}}} )$$, where I is the value of current in the ‘vertical’ frame, $${\delta ^{{{\hat{\Gamma }}_1}}}$$ and $${\delta ^{{{\hat{\Gamma }}_2}}}$$ are the delta-functions concentrated on the surfaces $${\hat{\Gamma }_1}$$ and $${\hat{\Gamma }_2}$$ (i.e. on the planes z = z1 and z = z2 corresponding to top and bottom frame borders, see Fig. 1) where the radial currents of the top and bottom NGR sources flow. Thus, the solution of both eqs (15) and (16) can be found as a result of solving one axisymmetrical problem: find the solution $$A_r^{}( {r,z} )$$ of equation   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta A_r^{} + \frac{1}{{{\mu _0}{r^2}}}A_r^{} + i\omega \sigma _r^{p,rz}A_r^{} = \frac{I}{{2\pi r}}( {{\delta ^{{{\hat{\Gamma }}_1}}} - {\delta ^{{{\hat{\Gamma }}_2}}}} ), \end{equation} (18)where $$\sigma _r^{p,rz}$$ is the lateral conductivity of a horizontally layered medium. The field $${{{\bf \vec{A}}}^{{\rm{I1,hor}}}}$$ defined by eq. (17) is not axisymmetrical. The vector function $${{{\bf \vec{J}}}^{{\rm{I1,}}{{\rm{P}}_1}{{\rm{P}}_2}}} - {{{\bf \vec{J}}}^{{\rm{I1,}}{{\rm{P}}_3}{{\rm{P}}_4}}}$$ defining the right-hand side of this equation has only one non-zero component depending on the orientation of the ‘vertical’ current frame: if the frame is oriented in the way shown in Fig. 1, then the y-component is the only non-zero component of the right-hand side vector of eq. (17) and, consequently, the vector function $${{{\bf \vec{A}}}^{{\rm{I1,hor}}}}$$ has the only non-zero component $$A_y^{{\rm{hor}}}( {x,y,z} )$$. If the ‘vertical’ current frame is oriented so that the lines P1P2 and P3P4 are parallel to the x-axis, then the x-component is the only non-zero component of the right-hand side vector of eq. (17) and, therefore, the vector function $${{{\bf \vec{A}}}^{{\rm{I1,hor}}}}$$ has the only non-zero component $$A_x^{{\rm{hor}}}( {x,y,z} )$$. These components can be found by solving the 3-D scalar equation   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta A_\xi ^{{\rm{hor}}} + i\omega \sigma _r^{p,rz}A_\xi ^{{\rm{hor}}} = I\left( {{\delta ^{{{\rm{P}}_1}{{\rm{P}}_2}}} - {\delta ^{{{\rm{P}}_3}{{\rm{P}}_4}}}} \right), \end{equation} (19)where $$A_\xi ^{{\rm{hor}}}( {x,y,z} )$$ is the vector potential component corresponding to the horizontal line direction (x-direction or y-direction), $${\delta ^{{{\rm{P}}_1}{{\rm{P}}_2}}}$$ and $${\delta ^{{{\rm{P}}_3}{{\rm{P}}_4}}}$$ are delta-functions concentrated on the lines P1P2 and P3P4. The right-hand side of eq. (19) can be approximated with any required accuracy by the set of point sources positioned along the lines P1P2 and P3P4 (see Fig. 1g) and, therefore, the 3-D component $$A_\xi ^{{\rm{hor}}}( {x,y,z} )$$ of the field induced by the sources concentrated on the lines P1P2 and P3P4 can be computed as a sum of axisymmetrical fields $$A_\xi ^{}( {r,z} )$$ induced by the point sources positioned along these lines. The field induced by two point sources (lying one below (above) another) can be calculated by the solution of the equation   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta A_\xi ^{} + i\omega \sigma _r^{p,rz}A_\xi ^{} = I{h_L}\left( {{\delta ^{{z_1}}} - {\delta ^{{z_2}}}} \right), \end{equation} (20)where $$A_\xi ^{}( {r,z} )$$ is the vector potential component corresponding to the ‘current line’ direction (i.e. x- or y-component), $${\delta ^{{z_1}}}$$ and $${\delta ^{{z_2}}}$$ are the delta-functions concentrated in the points (centres of the ‘current line’ subintervals, in which the point sources are located), hL is the length of the ‘current line’ subinterval (see Fig. 1g), and, I as before, is the value of the current in the frame. By summing axisymmetrical fields $$A_\xi ^{}( {r,z} )$$ taking into account the translation of the cylindrical coordinate system centres into the corresponding points of one of the ‘current lines’ of source ‘I1’, we obtain the 3-D distribution of one component of the vector potential, which describes the field of the ‘current line’. Thus, we can calculate the field induced by the non-axisymmetric source ‘I1’ by summing axisymmetrical fields taking into account the translation of the cylindrical coordinate system centres into corresponding points. The formulae for calculating the electric and magnetic field will be presented below (we take into account the fact that, when summing the components of the fields induced by separate parts of the source, the radial currents are compensated, and, therefore, the fields induced by these radial currents are not added to the resulting electric and magnetic field). In order to compute the axisymmetrical field induced by the source of type ‘I2’, we use the following mathematical model. Taking into account that the divergence of current density $${{{\bf \vec{J}}}^{{\rm{I2}}}} = {( {{J_r},{J_z}} )^{\rm{T}}}$$ is equal to zero, the electromagnetic field induced by each source ‘I2’ can be obtained by the solution of the differential equation system   \begin{equation} \! - \frac{{\rm{1}}}{{{\mu _0}}}\Delta A_r^{} + \frac{{\rm{1}}}{{{\mu _0}{r^2}}}A_r^{} + \sigma _r^{p,rz}\left( {i\omega A_r^{} + \frac{{\partial V}}{{\partial r}}} \right) = {J_r}, \end{equation} (21)  \begin{equation} \! - \frac{{\rm{1}}}{{{\mu _0}}}\Delta A_z^{} + \sigma _z^{p,rz}\left( {i\omega A_z^{} + \frac{{\partial V}}{{\partial z}}} \right) = {J_z}, \end{equation} (22)  \begin{equation} \! - \nabla \cdot \left( {{\sigma ^{p,rz}}\nabla V} \right) - i\omega \nabla \cdot \left( {{\sigma ^{p,rz}}{{{{\bf \vec{A}}}}^{{\rm{I2}}}}} \right) = 0, \end{equation} (23)where $${{{\bf \vec{A}}}^{{\rm{I2}}}} = {( {{A_r}( {r,z} ),{A_z}( {r,z} )} )^{\rm{T}}}$$ is the vector potential in the cylindrical coordinate system, V = V(r, z) is the scalar potential, and Δ = ∇ · ∇ is the scalar Laplacian. We represent the r-component of the vector potential $${{{\bf \vec{A}}}^{{\rm{I2}}}}$$ as $$A_r^{} = A_r^* + \tilde{A}_r^{}$$ where $$A_r^*$$ is the solution of equation in the form of (18). Then system (21)–(23) takes the form   \begin{equation} \! - \frac{{\rm{1}}}{{{\mu _0}}}\Delta \tilde{A}_r^{} + \frac{{\rm{1}}}{{{\mu _0}{r^2}}}\tilde{A}_r^{} + \sigma _r^{p,rz}\left( {i\omega \tilde{A}_r^{} + \frac{{\partial V}}{{\partial r}}} \right) = 0, \end{equation} (24)  \begin{equation} \! - \frac{{\rm{1}}}{{{\mu _0}}}\Delta A_z^{} + \sigma _z^{p,rz}\left( {i\omega A_z^{} + \frac{{\partial V}}{{\partial z}}} \right) = {J_z}, \end{equation} (25)  \begin{equation} \! - \nabla \cdot \left( {{\sigma ^{p,rz}}\nabla V} \right) - i\omega \nabla \cdot \left( {{\sigma ^{p,rz}}{{( {\tilde{A}_r^{},A_z^{}} )}^{\rm{T}}}} \right) = i\omega \nabla \cdot \left( {{\sigma ^{p,rz}}{{( {A_r^*,0} )}^{\rm{T}}}} \right). \end{equation} (26)The field $$A_r^*$$ in the right-hand side of eq. (26) is fully identical with $$A_r^{}$$ calculated for the source ‘I1’ (with the opposite sign) that allows us to solve eq. (18) only once and (what is more important!) to exclude these two mutually subtracting fields when we form the total field of the ‘vertical’ frame. Thus, in order to compute the field of the ‘vertical’ frame, it is necessary to solve three problems: eq. (18), (20) and equation system (24)–(26). Besides, the proposed scheme can be further simplified (that will increase its computational effectiveness) because there is no need to solve eq. (18). We define function u as $$u = \nabla \cdot {(A_r^*,0)^{\rm{T}}}$$. Differentiating eq. (18) properly we obtain the equation   \begin{equation} \! - \frac{1}{{{\mu _0}}}\Delta u + i\omega \sigma _r^{p,rz}u = I\left( {{\delta ^{P,1}} - {\delta ^{P,2}}} \right), \end{equation} (27)which is identical to eq. (20). Equation system (24)–(26) takes the form   \begin{equation} \! - \frac{{\rm{1}}}{{{\mu _0}}}\Delta \tilde{A}_r^{} + \frac{{\rm{1}}}{{{\mu _0}{r^2}}}\tilde{A}_r^{} + \sigma _r^{p,rz}\left( {i\omega \tilde{A}_r^{} + \frac{{\partial V}}{{\partial r}}} \right) = 0, \end{equation} (28)  \begin{equation} \! - \frac{{\rm{1}}}{{{\mu _0}}}\Delta A_z^{} + \sigma _z^{p,rz}\left( {i\omega A_z^{} + \frac{{\partial V}}{{\partial z}}} \right) = {J_z}, \end{equation} (29)  \begin{equation} \! - \nabla \cdot \left( {{\sigma ^{p,rz}}\nabla V} \right) - i\omega \nabla \cdot \left( {{\sigma ^{p,rz}}{{( {\tilde{A}_r^{},A_z^{}} )}^{\rm{T}}}} \right) = i\omega \sigma _r^{p,rz}u. \end{equation} (30)Finally, in order to obtain the field of the ‘vertical’ frame, it is necessary to solve two problems: eq. (20) and equation system (28)–(30). Further, we consider how to calculate the components of $${{{\bf \vec{E}}}^{p,2}}$$ using the solution of eq. (20) and equation system (28)–(30). In the case when the ‘vertical’ frame axis is parallel to the x-axis (i.e. the coordinates of its vertexes are (xc, y1, z1), (xc, y2, z1), (xc, y2, z2), and (xc, y1, z2)), the values of the components of $${{{\bf \vec{E}}}^{p,2}}$$ at the point (x, y, z) can be calculated with the help of the following formulae:   \begin{equation} E_x^{p,2}\left( {x,y,z} \right) = \left( {i\omega \tilde{A}_r^{}\left( {{r_1},z} \right) + \frac{{\partial V\left( {{r_1},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {x - {x_c}} \right)}}{{{r_1}}} - \left( {i\omega \tilde{A}_r^{}\left( {{r_2},z} \right) + \frac{{\partial V\left( {{r_2},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {x - {x_c}} \right)}}{{{r_2}}}, \end{equation} (31)  \begin{equation} E_y^{p,2}\left( {x,y,z} \right) = - i\omega A_y^{{\rm{hor}}} + \left( {i\omega \tilde{A}_r^{}\left( {{r_1},z} \right) + \frac{{\partial V\left( {{r_1},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {y - {y_1}} \right)}}{{{r_1}}} - \left( {i\omega \tilde{A}_r^{}\left( {{r_2},z} \right) + \frac{{\partial V\left( {{r_2},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {y - {y_2}} \right)}}{{{r_2}}}, \end{equation} (32)  \begin{equation} E_z^{p,2}\left( {x,y,z} \right) = \left( {i\omega \tilde{A}_z^{}\left( {{r_1},z} \right) + \frac{{\partial V\left( {{r_1},z} \right)}}{{\partial z}}} \right) - \left( {i\omega \tilde{A}_z^{}\left( {{r_2},z} \right) + \frac{{\partial V\left( {{r_2},z} \right)}}{{\partial z}}} \right), \end{equation} (33)where $${r_1} = \sqrt {{{( {x - {x_c}} )}^2} + {{( {y - {y_1}} )}^2}} $$, $${r_2} = \sqrt {{{( {x - {x_c}} )}^2} + {{( {y - {y_2}} )}^2}} $$ and $$A_y^{{\rm{hor}}}( {x,y,z} )$$ is obtained by solving eq. (27) as a result of summing the functions u(r, z) (after recalculating u(r, z) from the cylindrical coordinate system into the Cartesian one), which are calculated for the sources positioned in the centres of subintervals of the horizontal ‘current lines’ of the frame. The components of $${{{\bf \vec{B}}}^{p,2}}$$ at the point (x, y, z) can be calculated with the help of the following formulae:  \begin{equation*} B_x^{p,2}\left( {x,y,z} \right) = - \frac{{\partial\! A_y^{{\rm{hor}}}\left( {r,z} \right)}}{{\partial z}} + \left( {\frac{{\partial {A_r}\left( {{r_1},z} \right)}}{{\partial z}} - \frac{{\partial {A_z}\left( {{r_1},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {y - {y_1}} \right)}}{{{r_1}}} - \left( {\frac{{\partial {A_r}\left( {{r_2},z} \right)}}{{\partial z}} - \frac{{\partial {A_z}\left( {{r_2},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {y - {y_2}} \right)}}{{{r_2}}}, \end{equation*}   \begin{equation*} B_y^{p,2}\left( {x,y,z} \right) = - \left( {\frac{{\partial {A_r}\left( {{r_1},z} \right)}}{{\partial z}} - \frac{{\partial {A_z}\left( {{r_1},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {x - {x^c}} \right)}}{{{r_1}}} + \left( {\frac{{\partial {A_r}\left( {{r_2},z} \right)}}{{\partial z}} - \frac{{\partial {A_z}\left( {{r_2},z} \right)}}{{\partial r}}} \right) \cdot \frac{{\left( {x - {x^c}} \right)}}{{{r_2}}}, \end{equation*}   \begin{equation*} B_z^{p,2}\left( {x,y,z} \right) = \frac{{\partial\! A_y^{{\rm{hor}}}\left( {r,z} \right)}}{{\partial x}}. \end{equation*} If the transmitter coil is oriented so that the angle between its axis and the axis z is equal to α (see Fig. 2a), then $${{{\bf \vec{E}}}^p}$$ can be calculated using $${{{\bf \vec{E}}}^{p,1}}$$ and $${{{\bf \vec{E}}}^{p,2}}$$ by formula  \begin{equation*} {{{\bf \vec{E}}}^p} = {{{\bf \vec{E}}}^{p,1}}\cos \alpha + {{{\bf \vec{E}}}^{p,2}}{\sin \alpha}. \end{equation*} Figure 2. View largeDownload slide (a) The orientation of transmitter whose axis is perpendicular to the y-axis and (b) the rotation of the coordinate system when the transmitter orientation is arbitrary. Figure 2. View largeDownload slide (a) The orientation of transmitter whose axis is perpendicular to the y-axis and (b) the rotation of the coordinate system when the transmitter orientation is arbitrary. For arbitrary transmitter coil orientation, at first, we can calculate $${{{\bf \vec{\tilde{E}}}}^p}$$ by formulae (31)–(33) in the coordinate system $$\{ {\tilde{x},\tilde{y},\tilde{z}} \}$$, in which the coil axis is perpendicular to the $$\tilde{y}$$-axis, and then recalculate $${{{\bf \vec{\tilde{E}}}}^p}$$ into the coordinate system {x, y, z} by formulae $${{{\bf \vec{E}}}^p} = {{\bf M}}{{{\bf \vec{\tilde{E}}}}^p}$$ where M is the rotation matrix defined as follows:   \begin{equation} {{\bf M}} = \left( {\begin{array}{c@{\quad}c@{\quad}c} {\cos \beta }&{ - \sin \beta }&0\\ {\sin \beta }&{\cos \beta }&0\\ 0&0&1 \end{array}} \right). \end{equation} (34)In formula (34), β is the angle of rotation around the axis $$z = \tilde{z}$$ which transfers the coordinate system {x, y} to the coordinate system $$\{ {\tilde{x},\tilde{y}} \}$$(see Fig. 2b). The components of the magnetic field induction are calculated in a similar way. Thus, the proposed approach together with the possibility to take the field of the source in the horizontally layered medium as $${{{\bf \vec{E}}}^p}$$, which is much similar to the realistic 3-D medium in the vicinity of the logging tool location, makes it possible to exclude the appearance of the opposite sign components of the field, the moduli of which are close to each other. 3 FINITE-ELEMENT APPROXIMATION FE solutions of all the equations considered above and describing the secondary field $${{{\bf \vec{E}}}^s}$$ and the components of the primary field (which are axisymmetrical, if the coordinate system is properly translated) are found in the limited 3-D domain Ω (for $${{{\bf \vec{E}}}^s}$$) or 2-D domain Ωrz (for components of $${{{\bf \vec{E}}}^p}$$). On the outer boundary of Ω and Ωrz, we set homogeneous essential boundary conditions. For $${{{\bf \vec{E}}}^s}$$, these conditions represent the equality of $${{{\bf \vec{E}}}^s}$$ components to the zero, which are tangential to the boundary of Ω: $${{{\bf \vec{E}}}^s} \times {{\bf \vec{n}}} = {{\bf \vec{0}}}$$ where $${{\bf \vec{n}}}$$ is the normal to the boundary of Ω (these essential boundary conditions are set when the vector FEM is used). If the field $${{{\bf \vec{E}}}^s}$$ has a symmetry plane, then on this plane, we set either the homogeneous natural boundary condition $$( {\nabla \times {{{{\bf \vec{E}}}}^s}} ) \times {{\bf \vec{n}}} = {{\bf \vec{0}}}$$ (where $${{\bf \vec{n}}}$$ is the normal to the symmetry plane), if the current frame lies in this plane, or the homogeneous essential boundary condition $${{{\bf \vec{E}}}^s} \times {{\bf \vec{n}}} = {{\bf \vec{0}}}$$, if the current frame is sectioned by the symmetry plane. For calculating the axisymmetrical components of the field $${{{\bf \vec{E}}}^p}$$, on the outer boundaries of the computational domain Ωrz, we set homogeneous Dirichlet boundary conditions for all potentials describing the corresponding field. The size of computational domains Ω and Ωrz depends on the frequency and ‘average’ conductivity of the medium. Besides, the size of Ω is much smaller than the size of Ωrz. For example, for the frequency of 200 kHz, Ω = [−20,20] × [−15,15] × [−20,20] m3 (or Ω = [−20,20] × [0,15] × [−20,20] m3 if y = 0 is the symmetry plane of the field), and Ωrz = [0,70] × [−70,70] m2 (if the ‘average’ conductivity of the medium is about 0.05–1 Sm m−1). 3.1 Finite-element approximation for mathematical models used for $${{{\bf \vec{E}}}^p}$$ calculation The equivalent variational formulation for eq. (3), which defines the field of ‘horizontal’ coil (loop), takes the following form:   \begin{equation} \frac{1}{{{\mu _0}}}\left( {\int_{{{\Omega ^{rz}}}}{{\nabla E_\varphi ^{} \cdot \nabla \nu r{\rm{d}}r{\rm{d}}z}} + \int_{{{\Omega ^{rz}}}}{{\frac{1}{r}E_\varphi ^{}\nu {\rm{d}}r{\rm{d}}z}}} \right) + i\omega \int_{{{\Omega ^{rz}}}}{{\sigma _r^{p,rz}E_\varphi ^{}\nu r{\rm{d}}r{\rm{d}}z}} = - i\omega I\nu \left( {{r_\delta },{z_\delta }} \right){r_\delta }, \end{equation} (35)where ν(r, z) is a real test function, Ωrz is the computational domain in the cylindrical coordinate system. The expression in the right-hand side of eq. (35) is obtained as a result of the delta-function $${\delta ^{{J_\varphi }}}$$ (which concentrates on the circle r = rδ, z = zδ) operation on the test function ν(r, z). More specifically, the function $$E_\varphi ^{}$$ satisfying eqs (3) and (35) goes to infinity when approaching the circle r = rδ, z = zδ. The FE solution at the point (rδ, zδ) of the cylindrical domain Ωrz is finite, but together with mesh refinement around this point, the FE solution will grow infinitely. However, at some fixed (non-zero) distance from the point (rδ, zδ), the value of $$E_\varphi ^{}$$ can be calculated with any required accuracy by means of FEM. Therefore, when calculating $${{{\bf \vec{E}}}^{p,1}}$$, this fact must be taken into account, and when forming $${{{\bf \vec{E}}}^p}$$ for the solution of eq. (1), it is not allowed to come to the circle r = rδ, z = zδ closer than some pre-set distance. The solution $$E_\varphi ^{}$$ of eq. (35) is expressed as   \begin{equation} E_\varphi ^{} = \sum\limits_{j = 1}^{{n^{rz}}} {q_j^{{E_\varphi }}{\psi _j}} , \end{equation} (36)where ψj(r, z) are (real) nodal basis functions, nrz is their number (which coincides with the number of nodes in the 2-D mesh). As a result of substituting (36) into (35) and replacing the test function ν with all basis functions ψl in turn, the following system of linear algebraic equations (SLAE) is obtained:   \begin{equation} \!\left( {{{{\bf G}}^{{E_\varphi }}} + i\omega {{{\bf M}}^{{E_\varphi }}}} \right){{{\bf q}}^{{E_\varphi }}} = - i\omega {{{\bf f}}^{{E_\varphi }}}, \end{equation} (37)where the components of the matrices $${{{\bf G}}^{{E_\varphi }}}$$, $${{{\bf M}}^{{E_\varphi }}}$$ and vector $${{{\bf f}}^{{E_\varphi }}}$$ are calculated by formulae  \begin{equation*} G_{lj}^{{E_\varphi }} = \frac{1}{{{\mu _0}}}\left( {\int_{{{\Omega ^{rz}}}}{}\nabla {\psi _l} \cdot \nabla {\psi _j}r{\rm{d}}r{\rm{d}}z + \int_{{{\Omega ^{rz}}}}{}\frac{1}{r}{\psi _l}{\psi _j}{\rm{d}}r{\rm{d}}z} \right),\quad M_{lj}^{{E_\varphi }} = \int_{{{\Omega ^{rz}}}}{}\sigma _r^{p,rz}{\psi _l}{\psi _j}r{\rm{d}}r{\rm{d}}z, \end{equation*}   \begin{equation} f_l^{{E_\varphi }} = I{\psi _l}\left( {{r_\delta },{z_\delta }} \right){r_\delta }. \end{equation} (38)When FE approximation is used for this problem, we exclude a small subdomain in the vicinity of r = 0 from the computational domain Ωrz, on its boundary Γ0 the homogeneous Dirichlet boundary condition $${ {E_\varphi ^{}} |_{{\Gamma _0}}} = 0$$ is set (as on the distant boundaries of Ωrz). As mentioned above, for obtaining the field of the ‘vertical’ frame, it is necessary to solve two problems: eq. (20) and equation system (28)–(30). The equivalent variational formulation for eq. (20) takes the form   \begin{equation} \frac{1}{{{\mu _0}}}\int_{{{\Omega ^{rz}}}}{{\nabla A_\xi ^{} \cdot \nabla \nu r{\rm{d}}r{\rm{d}}z}} + i\omega \int_{{{\Omega ^{rz}}}}{{\sigma _r^{p,rz}A_\xi ^{}\nu r{\rm{d}}r{\rm{d}}z}} = Ih_L^{}\left( {\nu \left( {0,{z_1}} \right) - \nu \left( {0,{z_2}} \right)} \right), \end{equation} (39)where ν(0, z1) and ν(0, z2) are the results of delta-functions $${\delta ^{{z_1}}}$$ and $${\delta ^{{z_2}}}$$ operations on the test function ν(r, z). When this problem is solved in the vicinities of the points (0, z1) and (0, z2), FE solution can grow infinitely, and, therefore, when calculating $${{{\bf \vec{E}}}^{p,2}}$$, it is not allowed to come to these points closer than some pre-set distance. The solution $$A_\xi ^{}$$ of eq. (39) is also expressed as a linear combination of basis functions ψj(r, z)   \begin{equation} A_\xi ^{} = \sum\limits_{j = 1}^{{n^{rz}}} {q_j^{A_\xi ^{}}{\psi _j}} . \end{equation} (40)As a result of substituting (40) into (39) and replacing the test function with all basis functions ψl in turn, the following SLAE is obtained:   \begin{equation} \left( {{{{\bf G}}^{A_\xi ^{}}} + i\omega {{{\bf M}}^{A_\xi ^{}}}} \right){{{\bf q}}^{A_\xi ^{}}} = {{{\bf f}}^{A_\xi ^{}}}, \end{equation} (41)where the components of matrices $${{{\bf G}}^{A_\xi ^{}}}$$, $${{{\bf M}}^{A_\xi ^{}}}$$, and vector $${{{\bf f}}^{A_\xi ^{}}}$$ are calculated by formulae  \begin{equation*} G_{lj}^{A_\xi ^{}} = \frac{1}{{{\mu _0}}}\int_{{{\Omega ^{rz}}}}{}\nabla {\psi _l} \cdot \nabla {\psi _j}r{\rm{d}}r{\rm{d}}z,\quad M_{lj}^{A_\xi ^{}} = \int_{{{\Omega ^{rz}}}}{}\sigma _r^{p,rz}{\psi _l}{\psi _j}r{\rm{d}}r{\rm{d}}z, \end{equation*}   \begin{equation} f_l^{A_\xi ^{}} = Ih_L^{}\left( {{\psi _l}\left( {0,{z_1}} \right) - {\psi _l}\left( {0,{z_2}} \right)} \right). \end{equation} (42)The equivalent variational formulation for the equation system (28)–(30) takes the form   \begin{equation} \frac{1}{{{\mu _0}}}\left( {\int_{{{\Omega ^{rz}}}}{{\nabla \tilde{A}_r^{} \cdot \nabla \nu r{\rm{d}}r{\rm{d}}z}} + \int_{{{\Omega ^{rz}}}}{{\frac{1}{r}\tilde{A}_r^{}\nu {\rm{d}}r{\rm{d}}z}}} \right) + i\omega \int_{{{\Omega ^{rz}}}}{{\sigma _r^{p,rz}\tilde{A}_r^{}\nu r{\rm{d}}r{\rm{d}}z}} + \int_{{{\Omega ^{rz}}}}{{\sigma _r^{p,rz}\frac{{\partial V}}{{\partial r}}\nu r{\rm{d}}r{\rm{d}}z}} = 0, \end{equation} (43)  \begin{equation} \frac{1}{{{\mu _0}}}\int_{{{\Omega ^{rz}}}}{{\nabla A_z^{} \cdot \nabla \nu r{\rm{d}}r{\rm{d}}z}} + i\omega \int_{{{\Omega ^{rz}}}}{{\sigma _z^{p,rz}A_z^{}\nu r{\rm{d}}r{\rm{d}}z}} + \int_{{{\Omega ^{rz}}}}{{\sigma _z^{p,rz}\frac{{\partial V}}{{\partial z}}\nu r{\rm{d}}r{\rm{d}}z}} = \int_{{{\Omega ^{rz}}}}{{{J_z}\nu r{\rm{d}}r{\rm{d}}z}}, \end{equation} (44)  \begin{equation} \int_{{{\Omega ^{rz}}}}{{{\sigma ^{p,rz}}\nabla V \cdot \nabla \nu \,\,r{\rm{d}}r{\rm{d}}z}} + i\omega \int_{{{\Omega ^{rz}}}}{{{\sigma ^{p,rz}}{{( {\tilde{A}_r^{},A_z^{}} )}^{\rm{T}}} \cdot \nabla \nu \,\,r{\rm{d}}r{\rm{d}}z}} = i\omega \int_{{{\Omega ^{rz}}}}{{\sigma _r^{p,rz}u\nu r{\rm{d}}r{\rm{d}}z}}, \end{equation} (45)where ν(r, z) is the test function as before. Functions $$\tilde{A}_r^{}$$, $$A_z^{}$$ and V are represented in the form   \begin{equation} \tilde{A}_r^{} = \sum\limits_{i = 1}^{{n^{rz}}} {q_i^{\tilde{A}_r^{}}{\psi _i}} ,\quad A_z^{} = \sum\limits_{i = 1}^{{n^{rz}}} {q_i^{{A_z}}{\psi _i}} ,\quad V = \sum\limits_{i = 1}^{{n^{rz}}} {q_i^V{\psi _i}} . \end{equation} (46)Values $$q_i^{\tilde{A}_r^{}}$$, $$q_i^{{A_z}}$$ and $$q_i^V$$ are found from the SLAE which is obtained by substituting (46) into equation system (43)–(45) and replacing the test function ν with the basis function ψi in each equation. The obtained SLAE takes the following form:   \begin{equation} {{{\bf C}}^{AV}}{{{\bf q}}^{AV}} = {{{\bf f}}^{AV}}. \end{equation} (47)The size of system (47) is equal to 3nrz × 3nrz. The matrix CAV has a block structure with 3 × 3 block elements, vectors qAV and fAV also have a block structure with 3 × 1 block elements. The vector qAV has the following form: $${{{\bf q}}^{AV}} = {\{ {q_1^{\tilde{A}_r^{}},q_1^{A_z^{}},q_1^V,q_2^{\tilde{A}_r^{}},q_2^{A_z^{}},q_2^V,...} \}^{\rm{T}}}$$. The block elements of the matrix CAV take the following form:   \begin{equation} {{\bf C}}_{lj}^{AV} = \left[ {\begin{array}{c@{\quad}c@{\quad}c} {G_{ij}^{{{\tilde{A}}_r}} + i\omega M_{lj}^r}&0&{D_{lj}^r}\\ 0&{G_{ij}^{{A_z}} + i\omega M_{lj}^z}&{D_{lj}^z}\\ {i\omega D_{jl}^r}&{i\omega D_{jl}^z}&{G_{ij}^V} \end{array}} \right], \end{equation} (48)where   \begin{equation} G_{lj}^{{{\tilde{A}}_r}} = \frac{1}{{{\mu _0}}}\left( {\int_{{{\Omega ^{rz}}}}{}\nabla {\psi _l} \cdot \nabla {\psi _j}r{\rm{d}}r{\rm{d}}z + \int_{{{\Omega ^{rz}}}}{}\frac{1}{r}{\psi _l}{\psi _j}{\rm{d}}r{\rm{d}}z} \right),\quad G_{lj}^{{A_z}} = \frac{1}{{{\mu _0}}}\int_{{{\Omega ^{rz}}}}{}\nabla {\psi _l} \cdot \nabla {\psi _j}r{\rm{d}}r{\rm{d}}z, \end{equation} (49)  \begin{equation} G_{lj}^V = \int_{{{\Omega ^{rz}}}}{{{\sigma ^{p,rz}}\nabla {\psi _l} \cdot \nabla {\psi _j}\,\,r{\rm{d}}r{\rm{d}}z}},\quad D_{lj}^r = \int_{{{\Omega ^{rz}}}}{}\sigma _r^{p,rz}{\psi _l}\frac{{\partial {\psi _j}}}{{\partial r}}r{\rm{d}}r{\rm{d}}z,\quad D_{lj}^z = \int_{{{\Omega ^{rz}}}}{}\sigma _z^{p,rz}{\psi _l}\frac{{\partial {\psi _j}}}{{\partial z}}r{\rm{d}}r{\rm{d}}z, \end{equation} (50)  \begin{equation} M_{lj}^r = \int_{{{\Omega ^{rz}}}}{}\sigma _r^{p,rz}{\psi _l}{\psi _j}r{\rm{d}}r{\rm{d}}z,\quad M_{lj}^z = \int_{{{\Omega ^{rz}}}}{}\sigma _z^{p,rz}{\psi _l}{\psi _j}r{\rm{d}}r{\rm{d}}z. \end{equation} (51)The block elements of the vector fAV take the following form:   \begin{equation} {{\bf f}}_l^{AV} = \left[ {\begin{array}{@{}*{1}{c}@{}} 0\\ {f_l^z}\\ {i\omega f_l^u} \end{array}} \right],\quad f_l^z = \int_{{{\Omega ^{rz}}}}{{{J_z}{\psi _l}r{\rm{d}}r{\rm{d}}z}},\quad f_l^u = \int_{{{\Omega ^{rz}}}}{{\sigma _r^{p,rz}u{\psi _l}r{\rm{d}}r{\rm{d}}z}}. \end{equation} (52)Note that when calculating $$f_l^u$$ in (52), the function u is found as $$u = \sum_{j = 1}^{{n^{rz}}} {q_j^u{\psi _j}} $$ by solving the equation, which is analogous to eq. (20) (or 39), that is, with the use of eqs (41) and (42) when $$h_L^{} = 1$$. Therefore, the component $$f_l^u$$ in (52) can be calculated as $$f_l^u = \sum_j {M_{lj}^rq_j^u} $$ where $$M_{lj}^r$$ is defined by (51). There is a further point to be made here. In order to obtain the complex symmetrical SLAE, every third row of the system must be divided by iω. 3.2 Finite-element approximation for mathematical models used for calculating the 3-D field in anisotropic medium The equivalent variational formulation for eq. (1) takes the following form:   \begin{eqnarray} \int_{\Omega }{{\frac{1}{{{\mu _0}}}}}( {\nabla \times {{{{\bf \vec{E}}}}^s}} ) \cdot ( {\nabla \times {{\bf \vec{\boldsymbol\psi }}}} ){\rm{d}}\Omega + i\omega \int_{\Omega }{{\sigma {{{{\bf \vec{E}}}}^s}}} \cdot {{\bf \vec{\boldsymbol\psi }}}{\rm{d}}\Omega - {\omega ^2}\int_{\Omega }{{\varepsilon {{{{\bf \vec{E}}}}^s}}} \cdot {{\bf \vec{\boldsymbol\psi }}}{\rm{d}}\Omega &=& - i\omega \int_{\Omega }{{\left( {\sigma - {\sigma ^p}} \right){{{{\bf \vec{E}}}}^p}}} \cdot {{\bf \vec{\boldsymbol\psi }}}{\rm{d}}\Omega + {\omega ^2}\int_{\Omega }{{\varepsilon {{{{\bf \vec{E}}}}^p}}} \cdot {{\bf \vec{\boldsymbol\psi }}}{\rm{d}}\Omega \,\,, \end{eqnarray} (53)where Ω is the computational domain (and dΩ = dxdydz), the test vector function $${{\bf \vec{\boldsymbol\psi }}}$$ and the unknown vector function $${{{\bf \vec{E}}}^s}$$ are the members of the space Hcurl(Ω) of the vector functions $${{\bf \vec{\Phi }}}$$ with continuous tangential components, which satisfy conditions $$\int_{\Omega }{{{{\bf \vec{\Phi }}} \cdot {{\bf \vec{\Phi }}}{\rm{d}}\Omega }} < \infty $$ and $$\int_{\Omega }{{( {\nabla \times {{\bf \vec{\Phi }}}} ) \cdot ( {\nabla \times {{\bf \vec{\Phi }}}} ){\rm{d}}\Omega }} < \infty $$. The function $${{{\bf \vec{E}}}^s}$$ is expressed as   \begin{equation} {{{\bf \vec{E}}}^s} = \sum\limits_{j = 1}^n {q_j^s{{\vec{\boldsymbol\psi }}_j}} , \end{equation} (54)where $${\vec{\boldsymbol\psi }_j}$$ are the basis vector functions, n is their number (which coincides with the number of the edges in the 3-D mesh). As a result of substituting (54) into (53) and replacing the test function $${{\bf \vec{\boldsymbol\psi }}}$$ with the basis functions $${\vec{\boldsymbol\psi }_l}$$ in turn, we obtain the following SLAE:   \begin{equation} \left( {{{\bf G}} + i\omega {{{\bf M}}^\sigma } - {\omega ^2}{{{\bf M}}^\varepsilon }} \right){{{\bf q}}^s} = {{\bf f}}, \end{equation} (55)where the components of matrices G, Mσ and Mε, and the vector f are defined by relations   \begin{equation} {G_{lj}} = \frac{1}{{{\mu _0}}}\int_{\Omega }{}\left( {\nabla \times {{\vec{\boldsymbol\psi }}_l}} \right) \cdot \left( {\nabla \times {{\vec{\boldsymbol\psi }}_j}} \right){\rm{d}}\Omega ,\quad M_{lj}^\sigma = \int_{\Omega }{}\sigma {\vec{\boldsymbol\psi }_l} \cdot {\vec{\boldsymbol\psi }_j}{\rm{d}}\Omega ,\quad M_{lj}^\varepsilon = \int_{\Omega }{}\varepsilon {\vec{\boldsymbol\psi }_l} \cdot {\vec{\boldsymbol\psi }_j}{\rm{d}}\Omega , \end{equation} (56)  \begin{equation} {f_l} = - i\omega \int_{\Omega }{}\left( {\sigma - {\sigma ^p}} \right){{{\bf \vec{E}}}^p} \cdot {\vec{\boldsymbol\psi }_l}{\rm{d}}\Omega + {\omega ^2}\int_{\Omega }{{\varepsilon {{{{\bf \vec{E}}}}^p}}} \cdot {\vec{\boldsymbol\psi }_l}{\rm{d}}\Omega . \end{equation} (57)The function $${{{\bf \vec{E}}}^p}$$ is expressed as a linear combination of the basis functions $${\vec{\boldsymbol\psi }_j}$$  \begin{equation} {{{\bf \vec{E}}}^p} = \sum\limits_{j = 1}^n {q_j^p{{\vec{\boldsymbol\psi }}_j}} , \end{equation} (58)where expansion coefficients $$q_j^p$$ in basis $$\{ {{{\vec{\boldsymbol\psi }}_j}} \}$$ are the values of tangential components of $${{{\bf \vec{E}}}^p}$$ on the FEs’ edges. In this case, the vector f in the right-hand side of FE equation system (55) can be calculated by the formula   \begin{equation} {{\bf f}} = \left( { - i\omega {{{{\bf \bar{M}}}}^\sigma } + {\omega ^2}{{{\bf M}}^\varepsilon }} \right){{\bf q}}_{}^p, \end{equation} (59)where components of the matrix $${{{\bf \bar{M}}}^\sigma }$$ are defined by the formula   \begin{equation} \bar{M}_{lj}^\sigma = \int_{\Omega }{}\left( {\sigma - {\sigma ^p}} \right){\vec{\boldsymbol\psi }_l} \cdot {\vec{\boldsymbol\psi }_j}{\rm{d}}\Omega . \end{equation} (60) 4 SPECIAL ASPECTS OF IMPLEMENTATION For the solution of axisymmetrical problems (when we calculate the primary field $${{{\bf \vec{E}}}^p}$$), we use the regular rectangular meshes with smoothly increasing of the step when moving away from the source, which is located on the z-axis. The step in the vicinity of the source is set about 0.01 m, the factor χ of increase of the mesh step when moving away from the source in the coordinates R and Z is set in the range of 1.05–1.08 (i.e. the value of the next mesh step hnext in each of the coordinates is defined by the value of the previous step hprev in the corresponding coordinate as hnext = χ ⋅ hprev). The approach proposed in the previous sections for modelling the primary field requires much lower computational costs in comparison with computational costs required for calculating the 3-D secondary field. Therefore, we use the special algorithms of optimization only for the 3-D meshes. These algorithms make it possible to reduce the computational time required for calculating the secondary field by approximately an order in comparison with the calculations on regular meshes (which are generated according to the principle considered above for 2-D meshes). Therefore, in this Section, we consider the algorithms of generating the optimized irregular non-conforming 3-D meshes and principles of constructing the FE approximation on them, as well as another aspects (grouping of logging tool positions and automatic choice of the background medium), which allow additional significant reduction of the computational costs. Note that we do not use the approaches with adaptive mesh generation, which allow the convergence towards the solution using the sequence of meshes. These approaches are based on various error estimators (Schwarzbach et al.2011; Grayver & Kolev 2015; Yin et al.2016, etc.) and require the additional costs caused by the calculations on several meshes in order to obtain the adaptive mesh and for each transmitter position. At the same time, these approaches do not always guarantee the required accuracy when the signals are calculated in the receivers. This fact is implicitly confirmed by Grayver & Kolev (2015) where, after generating the adaptive mesh (and satisfying the corresponding criteria), the additional refinements are performed in the vicinities of the receivers. 4.1 Generation of the optimized finite-element meshes For the FE solution, we use the optimized meshes with local refinements. Note that the trajectory length in the logging problems is usually much bigger than the size of the computational domain, in which the field for one position of the source can be computed with the required accuracy. For increasing the numerical solution accuracy it is necessary to refine the mesh in the vicinity of the transmitters and receivers. Therefore, during the calculation of the signals received by the logging tool when it is moving along the trajectory, it is necessary to change not only the positions of the local refinements of the mesh but also the computational domain itself (the analogous approach has been considered in (Yang et al.2014), but it uses regular meshes). The accuracy of approximation of the computational domain geometry is the important aspect for obtaining the high accuracy solution. On the one hand, this approximation should be quite good in the vicinity of the logging tool (i.e. the error caused by the inaccurate boundaries approximation should not exceed the error of the numerical solution approximation). On the other hand, the excessive mesh refinement caused by good approximation of the boundaries, which are distant from the logging tool, can lead to the significant increase in the number of FEs and, consequently, computational costs. The numerical results have shown that high accuracy geometry approximation is required only in some subdomain around the logging tool, and rough geometry approximation can be used at some distance from it. Note that though we use the irregular non-conforming meshes, we do not use FEs with curvilinear or inclined boundaries. This is related to the fact that, if we use the FEs with curvilinear or inclined boundaries, which smoothly approximate the corresponding boundaries of the domain, the application of the primary-secondary field approach can lead to the intersection of FEs of the 3-D mesh by the layers interfaces of the background medium. In this case, the secondary field $${{{\bf \vec{E}}}^s}$$ should be discontinuous inside those FEs which are intersected by the interfaces of background medium layers with different conductivity (because $${{{\bf \vec{E}}}^t} = {{{\bf \vec{E}}}^p} + {{{\bf \vec{E}}}^s}$$ inside the FE should be continuous, but $${{{\bf \vec{E}}}^p}$$ will be discontinuous if the FE is intersected by the layers with a different conductivity). Therefore, we use another approach (stair-like approximation) in which the FEs are never intersected by the layers interfaces of the background medium. Mesh generation algorithm consists of four main stages. At the first stage, the size and structure of the computational domain are defined depending on the logging tool position (or several positions; it will be described in Section 4.3). The second stage is the stage of generating the so-called base mesh, which does not include inner borders of the computational domain and is refined in the subdomain containing the logging tool. At the third stage, the geoelectrical inhomogeneities are embedded into the mesh taking into account the requirements to the accuracy of their boundaries approximation depending on the distance from the logging tool. And, finally, at the fourth stage, the local refinements in the vicinity of receivers and transmitters are performed. Thus, when generating the base mesh, we select the parallelepiped subdomain, which contains the transmitters and receivers locations for one tool position (or several tool positions included in one group). Inside this subdomain, the base mesh is generated with a constant step in all directions, and outside this subdomain, the cells size increases smoothly with moving away from the selected subdomain boundaries. Thus, the base mesh for a group of tool positions on the trajectory is defined by three parameters. They are the following: the step inside the selected subdomain, the distance from the tool to boundaries of the computational domain, and the factor of increase of the cell size (when moving away from the selected subdomain boundaries). Based on multiple numerical experiments performed for different geoelectrical models, we have obtained the following values of these parameters. For the frequency of 100 kHz, which we take as a ‘basic’ one, the size of the minimal step in the base mesh is set equal to 0.5 m, and the distance from the tool to the boundaries of the computational domain is set equal to 20 m in all directions. For other frequencies, these sizes are scaled by multiplying by the square root of the frequencies ratio. The factor of increase of the cells size (when moving away from the boundaries of the subdomain containing all locations of transmitters and receivers included in a group) is taken equal to 1.2. Before the third stage of the algorithm, all geoelectrical inhomogeneities included in the computational domain (or the whole computational domain) are approximated by parallelepipeds the size of which is defined by the step of geometry approximation. The step defining the approximation of geometry of embedded geoelectrical inhomogeneities in the vicinity of receivers’ locations is about 0.2 m for the problems without a borehole and invasion zone, and it is about 0.05 m for the problems in which the borehole and invasion zone need to be taken into account. Outside the vicinity of receivers’ locations, the step of boundaries approximation increases proportionally to the increase in the step size of the base mesh. The local refinements in the vicinity of receivers (at the fourth stage of the algorithm) are defined by the following parameters: the number of refinements inside the FEs into which the receivers get, the number of neighbouring levels of FEs, and the number of refinements for each neighbouring levels. For the examples considered in this paper, the values of these parameters are taken equal to 2, 3 and 1, respectively. As mentioned above, for solving the problems considered here, we use the mathematical models based on the representation of the total field as a sum of the primary (corresponding to the field induced by the source in the horizontally layered background medium) and secondary (which is the impact field of 3-D and 2-D geoelectrical inhomogeneities) components. Note that the local refinements in the vicinity of the transmitter should be performed only in the case if the transmitter is located in the subdomain (of the computational domain) in which σp ≠ σ. As an example, consider the geoelectrical model the vertical section of which is shown in Fig. 3. The logging tool trajectory is shown by black solid line. Figure 3. View largeDownload slide The geoelectrical model specific for induction logging problems. The model contains five layers (interfaces of which are slightly sloped) and the fault. The layers of equivalent conductivity are shown by the same colour. The conductivity of the middle layer is anisotropic: lateral conductivity is 0.256 Sm m−1 and vertical conductivity is 0.025 Sm m−1. Figure 3. View largeDownload slide The geoelectrical model specific for induction logging problems. The model contains five layers (interfaces of which are slightly sloped) and the fault. The layers of equivalent conductivity are shown by the same colour. The conductivity of the middle layer is anisotropic: lateral conductivity is 0.256 Sm m−1 and vertical conductivity is 0.025 Sm m−1. We consider the logging tool that includes one transmitter and two receivers the offsets of which are 4.5 and 9.5 m. The distance between two adjacent logging tool positions is 1 m. The value of the frequency is equal to 100 kHz. Fig. 4(a) shows the base mesh, which is obtained after the first and second stages of the algorithm, for one group of three neighbouring logging tool positions on the trajectory. Fig. 4(b) shows the mesh obtained after the third stage of the algorithm, that is, after embedding 2-D geoelectrical inhomogeneities into the mesh. Fig. 5(a) shows the mesh obtained after the fourth stage of the algorithm, that is, after local refinement in the vicinity of the receivers and transmitters. Fig. 5(b) shows the 3-D view of the sectioned FE mesh. In Fig. 6, the FE mesh is shown for other logging tool positions on the trajectory. Figure 4. View largeDownload slide (a) The vertical section of the base mesh obtained after the first and second algorithm stages and (b) the vertical section of the mesh obtained after embedding 2-D geoelectrical inhomogeneities into the mesh. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the ith transmitter. Figure 4. View largeDownload slide (a) The vertical section of the base mesh obtained after the first and second algorithm stages and (b) the vertical section of the mesh obtained after embedding 2-D geoelectrical inhomogeneities into the mesh. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the ith transmitter. Figure 5. View largeDownload slide (a) The vertical section of the mesh obtained after local refinements in the vicinities of receivers and transmitters and (b) 3-D view of the sectioned finite-element mesh. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the ith transmitter. Figure 5. View largeDownload slide (a) The vertical section of the mesh obtained after local refinements in the vicinities of receivers and transmitters and (b) 3-D view of the sectioned finite-element mesh. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the ith transmitter. Figure 6. View largeDownload slide The vertical sections of the meshes generated for different groups of the adjacent tool positions. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the jth transmitter. Figure 6. View largeDownload slide The vertical sections of the meshes generated for different groups of the adjacent tool positions. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the jth transmitter. 4.2 Finite-element approximation on the meshes with hanging nodes 4.2.1 ‘Non-conforming’ and ‘conforming’ bases on the meshes with hanging nodes The meshes presented in the previous section contain the non-conforming FEs. On these meshes, we introduce two sets of the basis functions: $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$ and $$\{ {\vec{\boldsymbol\psi }_k^c} \}$$, where the superscript ‘nc’ denotes non-conforming, and superscript ‘c’ denotes conforming. The basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$ are associated with all edges of the mesh with hanging nodes. On each FE Ωe containing the corresponding (kth) edge, $$\vec{\boldsymbol\psi }_k^{nc}$$ is defined in the standard way: at this edge, the basis function $$\vec{\boldsymbol\psi }_k^{nc}$$ has the only non-zero (equal to 1) component, which corresponds to the coordinate (x, y or z) coinciding with the edge direction, and, on Ωe, this component decreases to zero as a bilinear function of two other coordinates when coming to three other edges with the same direction. Obviously, the functions $$\vec{\boldsymbol\psi }_k^{nc} \notin {{{\bf H}}^{{\rm{curl}}}}( \Omega )$$ occur in the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. As an example, consider the fragment of the irregular non-conforming mesh (with hanging nodes), which is shown in Fig. 7. Denote the basis function, which is associated with the edge connecting nodes with numbers 5 and 26, as $$\vec{\boldsymbol\psi }_5^{nc}$$. This function is non-zero on the FEs Ω1, Ω2, and Ω3, and zero on all other FEs. Further, we will denote the edges, FEs, and their faces by the sets of nodes, which are listed in braces. Thus, the basis function $$\vec{\boldsymbol\psi }_5^{nc}$$ is associated with the edge {5, 26}, FE Ω1 = {1, 2, 4, 5, 22, 23, 25, 26}, FE Ω1 has the faces {1, 2, 4, 5}, {1, 2, 22, 23}, etc. As an example, consider the basis functions $$\vec{\boldsymbol\psi }_{36}^{nc}$$, $$\vec{\boldsymbol\psi }_{37}^{nc}$$ and $$\vec{\boldsymbol\psi }_{38}^{nc}$$, which are associated with the edges {5, 10}, {10, 16} and {16, 26}, respectively. These edges are the sub-edges of the edge {5, 26}. The basis function $$\vec{\boldsymbol\psi }_{36}^{nc}$$ is non-zero only on FE Ω4 = {5, 6, 8, 9, 10, 12, 13, 15}, basis function $$\vec{\boldsymbol\psi }_{37}^{nc}$$ is non-zero only on FE Ω5 = {10, 11, 13, 14, 16, 17, 19, 20}, basis function $$\vec{\boldsymbol\psi }_{38}^{nc}$$ is non-zero only Ω7 = {16, 17, 19, 20, 26, 27, 30, 31}, and on all other FEs, these basis functions are equal to zero. Figure 7. View largeDownload slide The fragment of the irregular non-conforming finite-element mesh with hanging nodes. Figure 7. View largeDownload slide The fragment of the irregular non-conforming finite-element mesh with hanging nodes. Besides the basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$ associated with the sub-edges of ‘big’ edges (which are the edges of other FEs), in the ‘non-conforming’ basis, there are functions, which are associated with the FEs’ edges lying within the faces of other FEs (naturally, such edges do not coincide with the boundaries of these faces). In Fig. 7, these are the functions $$\vec{\boldsymbol\psi }_{39}^{nc}$$ and $$\vec{\boldsymbol\psi }_{40}^{nc}$$, which are associated with the edges {11, 17} and {17, 27}, respectively. The function $$\vec{\boldsymbol\psi }_{39}^{nc}$$ is non-zero only on FE Ω5 (see above its nodes) and Ω6 = {11, 12, 14, 15, 17, 18, 20, 21}, and the function $$\vec{\boldsymbol\psi }_{40}^{nc}$$ is non-zero only on FE Ω7 and Ω8 = {17, 18, 20, 21, 27, 28, 31, 32}. In general, by the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$ we mean the basis functions associated with all edges (including sub-edges). These basis functions are defined in the standard form on the FEs including these edges. On the non-conforming FE mesh with hanging nodes, we introduce the basis functions $$\vec{\boldsymbol\psi }_j^c$$ with continuous tangential components in the computational domain Ω, that is, $$\vec{\boldsymbol\psi }_j^c \in {{{\bf H}}^{{\rm{curl}}}}( \Omega )$$. These basis functions are associated only with the ‘main’ edges of the mesh. By the ‘main’ edges we mean such edges that are not the sub-edges and do not lie within the faces of any FEs. Define the basis functions $$\vec{\boldsymbol\psi }_j^c$$ in the form of linear combination of the basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$  \begin{equation} \vec{\boldsymbol\psi }_j^c = \sum\limits_k {{T_{jk}}\vec{\boldsymbol\psi }_k^{nc}} , \end{equation} (61)that is, the basis $$\{ {\vec{\boldsymbol\psi }_j^c} \}$$ is constructed from the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$ with the help of the transition matrix T = {Tjk}. The size of this matrix is nc × n, where n is the total number of the edges in the mesh (including sub-edges and edges, which lie within the faces of any FEs), and nc is the number of ‘main’ edges only. Thus, specifically with the help of the transition matrix T, we construct the basis $$\{ {\vec{\boldsymbol\psi }_j^c} \}$$ in which all vector functions have continuous tangential components in Ω (i.e. ∀j$$\vec{\boldsymbol\psi }_j^c \in {{{\bf H}}^{{\rm{curl}}}}( \Omega )$$). We note a very important aspect related to the ‘non-conforming’ basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. Besides the basis functions with the discontinuous tangential components (which we have called ‘non-conforming’), the ‘conforming’ basis functions (with continuous tangential components) are also included into the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. They are associated with the edges, which we will name ‘regular’. These edges do not contain subedges, are not subedges themselves, and do not lie within any face. For example, in Fig. 7, the ‘regular’ edges are the edges {1, 22}, {4, 25}, {2, 23}, etc. (if we assume that either the faces of FE Ω1 lie on the computational domain boundary, or FE Ω1 and FEs adjacent to Ω1, which are not shown in this mesh fragment, are conforming). The superscript ‘nc’ added to these basis functions (associated with ‘regular’ edges) denotes only that they are included in the basis in which there are vector functions with discontinuous tangential components. Thus, the superscript ‘nc’ added to the function from the ‘non-conforming’ basis does not always denote that this function necessarily have discontinuous tangential components, that is, among the functions from $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$, there are functions with continuous tangential components, which are associated with the ‘regular’ edges. Therefore, further, the superscript ‘nc’ added to the basis function $$\vec{\boldsymbol\psi }_k^{nc}$$ will not rigorously denote that this function necessarily has discontinuous tangential components (i.e. the term ‘conforming’ function $$\vec{\boldsymbol\psi }_k^{nc}$$ is applicable and denotes that such a function from the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$ is associated with the ‘regular’ edge). However, all ‘regular’ edges are the ‘main’, but unlike ‘regular’, the ‘main’ edges can contain the sub-edges (which are the edges of some FEs). For example, in Fig. 7, the edge {5, 26} is ‘main’, but not ‘regular’, because it contains subedges {5, 10} (the edge of FE Ω4), {10, 16} (the edge of FE Ω5), and {16, 26} (the edge of FE Ω7). 4.2.2 The construction of the ‘conforming’ basis from the ‘non-conforming’ Consider how the ‘conforming’ basis function $$\vec{\boldsymbol\psi }_j^c$$ can be constructed from the ‘non-conforming’ functions $$\vec{\boldsymbol\psi }_k^{nc}$$. In so doing, we really develop the algorithm of generating the matrix T. Let the face of one FE (this face will be called ‘big’) be adjacent to the several faces of other FEs (these faces will be called ‘small’). Then, in order to ensure on this ‘big’ face the continuity of the tangential component of the basis function $$\vec{\boldsymbol\psi }_j^c$$ associated with one of its edges, it is necessary to define the weights Tjk in the expansion (61) in basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$. The weights Tjk are equal to the values of the tangential component $$\vec{\boldsymbol\psi }_j^c$$ on the corresponding edges of the ‘small’ faces. We illustrate this by the example of the mesh shown in Fig. 7. Consider the face {5, 6, 26, 28} of FE Ω2 (this is the ‘big’ face) adjacent to the faces {5, 6, 10, 12} of FE Ω4, {10, 11, 16, 17} of FE Ω5, {11, 12, 17, 18} of FE Ω6, {16, 17, 26, 27} of FE Ω7 and {17, 18, 27, 28} of FE Ω8 (these are the ‘small’ faces). Let, as before, the vector function associated with the edge {5, 26} have number 5 in the ‘non-conforming’ basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$, and vector functions associated with the edges {5, 10}, {10, 16}, {16, 26}, {11, 17} and {17, 27} have numbers 36, 37, 38, 39 and 40, respectively. In the ‘conforming’ basis $$\{ {\vec{\boldsymbol\psi }_j^c} \}$$, the vector function associated with the edge {5, 26} has number 5 too. Then, the ‘conforming’ vector function associated with the edge {5, 26} takes the form   \begin{equation} \vec{\boldsymbol\psi }_5^c = \vec{\boldsymbol\psi }_5^{nc} + \vec{\boldsymbol\psi }_{36}^{nc} + \vec{\boldsymbol\psi }_{37}^{nc} + \vec{\boldsymbol\psi }_{38}^{nc} + \frac{{{h_{11,12}}}}{{{h_{10,12}}}}\left( {\vec{\boldsymbol\psi }_{39}^{nc} + \vec{\boldsymbol\psi }_{40}^{nc}} \right), \end{equation} (62)where hν, η is the length of the edge {ν,η}. Coefficients T5, 5 = 1, T5, 36 = 1, T5, 37 = 1, T5, 38 = 1, $${T_{5,39}} = {{{h_{11,12}}} / {{h_{10,12}}}}$$, and $${T_{5,40}} = {{{h_{11,12}}} / {{h_{10,12}}}}$$ are the values of the function $$\vec{\boldsymbol\psi }_5^c$$ at the edges {5, 10}, {10, 16}, {16, 26}, {11, 17} and {17, 27}, which lie on the ‘big’ face {5, 6, 26, 28} and with which the functions $$\vec{\boldsymbol\psi }_{36}^{nc}$$, $$\vec{\boldsymbol\psi }_{37}^{nc}$$, $$\vec{\boldsymbol\psi }_{38}^{nc}$$, $$\vec{\boldsymbol\psi }_{39}^{nc}$$ and $$\vec{\boldsymbol\psi }_{40}^{nc}$$, are associated. Note that the function $$\vec{\boldsymbol\psi }_5^c$$ identically coincides with the function $$\vec{\boldsymbol\psi }_5^{nc}$$ on the ‘big’ face, because on FE Ω2, both $$\vec{\boldsymbol\psi }_5^c$$ and $$\vec{\boldsymbol\psi }_5^{nc}$$ are defined by the same shape function. We also note that formula (62) defines the function $$\vec{\boldsymbol\psi }_5^c$$ only inside the FEs, and its values on the faces are defined as continuous extensions of its values from the FEs. On the face {5, 6, 26, 28} itself, formula (62) should be interpreted only in the sense of continuous extensions of the basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$ onto this face from the FEs, which contain the edges corresponding to $$\vec{\boldsymbol\psi }_k^{nc}$$, because all functions $$\vec{\boldsymbol\psi }_k^{nc}$$ in the right-hand side of (62) are discontinuous on this face. As a result, we really obtain the function $$\vec{\boldsymbol\psi }_5^c$$ with continuous tangential components on the face {5, 6, 26, 28}. Thus, we propose the following algorithm of calculating the matrix T, which makes it possible to obtain the ‘conforming’ basis functions $$\vec{\boldsymbol\psi }_j^c$$ from the ‘non-conforming’ (but standard for the vector FEM) basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$. To begin with, we assume that in the ‘non-conforming’ basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$, the functions associated with the ‘main’ edges are numbered first. The same numbering of the ‘main’ edges is used for the functions from the ‘conforming’ basis $$\{ {\vec{\boldsymbol\psi }_j^c} \}$$. Then, the left block of the matrix T with the size nc × nc is the identity matrix (i.e. Tjj = 1, and Tjk = 0 if k ≠ j). This is convenient both for the storage of the matrix T (i.e. in fact, there is no need to store this block) and for assembling the matrices of FEs. Note that the basis functions $$\vec{\boldsymbol\psi }_j^c$$ associated with the ‘regular’ edges are identically equal to $$\vec{\boldsymbol\psi }_j^{nc}$$, and for these edges, all components of the corresponding (jth) rows in the right block of the matrix T are equal to zero (i.e. Tjk = 0 if k > nc). Non-zero components Tjk if k > nc are calculated in the following way. If the edge numbered k is the sub-edge of the edge numbered α ≤ nc, then Tαk = 1, and all other components of the kth column of the matrix T are zero (i.e. ∀β ≠ α Tβk = 0). If the edge numbered k lies within the ‘big’ face, on which the basis functions $$\vec{\boldsymbol\psi }_\alpha ^{nc}$$ and $$\vec{\boldsymbol\psi }_\beta ^{nc}$$ for α, β ≤ nc have non-zero component tangential to this edge, then only two components Tαk and Tβk are non-zero in the kth column of the matrix T, and their values are equal to the values of functions $$\vec{\boldsymbol\psi }_\alpha ^{nc}$$ and $$\vec{\boldsymbol\psi }_\beta ^{nc}$$ at the kth edge, respectively. Thus, for the mesh shown in Fig. 7, if the edge {11, 17}, which lies within the face {5, 6, 26, 28}, has number 39, and edges {5, 26} and {6, 28} have numbers 5 and 6, respectively, then only the components $${T_{5,39}} = {{{h_{11,12}}} / {{h_{10,12}}}}$$ and $${T_{6,39}} = {{{h_{10,11}}} / {{h_{10,12}}}}$$ (where, as before, hν, η is the length of the edge {ν,η}) are non-zero in the 39th column of the matrix T. Finally, we consider the situation when the edge numbered k > nc lies within the ‘big’ face and even if one of the edges of this ‘big’ face, which are parallel to the k-th edge, is not the ‘main’, that is, its number is α > nc. This situation is shown in Fig. 8. The edge {10, 16} lies within the face {4, 5, 24, 26} and, therefore, this face is ‘big’ over the edge {10, 16} and FEs {9, 10, 12, 13, 15, 16, 18, 19} and {10, 11, 13, 14, 16, 17, 19, 20}, which contain it. But the edge {4, 24} of this ‘big’ face lies within the face {2, 7, 22, 28}, which is ‘big’ over this edge and FE {2, 3, 4, 5, 22, 23, 24, 26}, which contains this edge. Therefore, the edge {4, 24} is not ‘main’ and its number is α > nc. In this case, the non-zero components of the kth column of the matrix T are calculated with the use of a ‘chain’. First, the values $${\tilde{T}_{\alpha k}}$$ are calculated as values of tangential components of the basis function $$\vec{\boldsymbol\psi }_\alpha ^{nc}$$ at the kth edge. For the example shown in Fig. 8 and edges {10, 16} (its number is equal to k) and {4, 24} (its number is equal to α), $${\tilde{T}_{\alpha k}}$$ is equal to $${{{h_{10,11}}} / {{h_{9,11}}}}$$. The tilde ‘∼’ over Tαk have been added because if α > nc, there are no such components in the matrix T (because its size is nc × n), and the values $${\tilde{T}_{\alpha k}}$$ if α > nc are really fictitious and used only in the ‘chain’ for calculating the non-fictitious components of the matrix T (i.e. components Tjk if j ≤ nc). Figure 8. View largeDownload slide The fragment of the irregular non-conforming finite-element mesh with hanging nodes such that the ‘chains’ appear when the transition matrix is generated. Figure 8. View largeDownload slide The fragment of the irregular non-conforming finite-element mesh with hanging nodes such that the ‘chains’ appear when the transition matrix is generated. Further, the ‘big’ face within which the edge numbered α lies is processed. For the edges of this face, which are parallel to the edge numbered α, the values of basis functions $$\vec{\boldsymbol\psi }_\gamma ^{nc}$$ and $$\vec{\boldsymbol\psi }_\nu ^{nc}$$ associated with these edges are calculated at the edge numbered α. If γ ≤ nc and ν ≤ nc, that is, the edges numbered γ and ν are the ‘main’ ones, then the components $${T_{\gamma k}} = {T_{\gamma \alpha }} \cdot {\tilde{T}_{\alpha k}}$$ and $${T_{\nu k}} = {T_{\nu \alpha }} \cdot {\tilde{T}_{\alpha k}}$$ are non-zero components in the kth column of the matrix T. Thus, for the mesh shown in Fig. 8, if the edges {10, 16}, {2, 22} and {7, 28} are numbered k, γ and ν, respectively, then $${T_{\gamma ,k}} = \frac{{{h_{10,11}}}}{{{h_{9,11}}}} \cdot \frac{{{h_{4,7}}}}{{{h_{2,7}}}}$$ and $${T_{\nu ,k}} = \frac{{{h_{10,11}}}}{{{h_{9,11}}}} \cdot \frac{{{h_{2,4}}}}{{{h_{2,7}}}}$$. In so doing, the generation of the kth column of the matrix T is completed. If at least one of the edges is not the ‘main’ one (i.e. γ > nc and/or ν > nc), then for this edge (or these two edges), instead of Tγα and/or Tνα the ‘fictitious’ values $${\tilde{T}_{\gamma \alpha }}$$ and/or $${\tilde{T}_{\nu \alpha }}$$ are calculated, and the corresponding ‘chains’ (or one of them) must be continued in a similar manner. Thus, the generation of the kth column of the matrix T will be completed when all ‘chains’, which have already been begun, are finished (i.e. the ‘main’ edge should be at the end of each ‘chain’). Note that a number of non-zero components of the matrix T (even without taking into account the identity submatrix with the size of nc × nc) is small. Therefore, for its storage it is reasonable to use the sparse format as for the global matrix of the FE equations system. The use of the matrix T allows sufficiently convenient generation of matrices G, Mσ, and Mε of FE equations system (55) for the ‘conforming’ basis $$\{ {\vec{\boldsymbol\psi }_j^c} \}$$ by assembling the standard FE matrices corresponding to the ‘non-conforming’ basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. 4.2.3 Generation of global matrices of finite-element equations system We find the solution $${{{\bf \vec{E}}}^s}$$ of eq. (53) as a linear combination of basis functions $$\vec{\boldsymbol\psi }_j^c$$:  \begin{equation*} {{{\bf \vec{E}}}^s} = \sum\limits_{j = 1}^{{n_c}} {q_j^s\vec{\boldsymbol\psi }_j^c} . \end{equation*} Note once again that basis functions $$\vec{\boldsymbol\psi }_j^c$$ have continuous tangential components in the whole computational domain Ω ($$\vec{\boldsymbol\psi }_j^c \in {{{\bf H}}^{{\rm{curl}}}}( \Omega )$$) and, therefore, are suitable for variational equation (53). We show how the matrices G, Mσ, and Mε, which define the matrix of FE equations system (55), can be calculated using the FE matrices, the components of which are obtained with the use of the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. The basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$ is defined by the standard shape vector functions on FEs, including the ‘non-conforming’ ones, that is, on each FE the basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$, which are non-zero on it, are the shape functions. Thus, shape functions are the local basis functions (on FE) of the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. We obtain the expressions for calculating the components of the matrix G (see eq. 56) with the use of the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. Considering that tangential components of basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$ can be discontinuous on the faces of the FEs Ωe, before turning to the basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$, we represent the integral over $$\Omega = \bigcup\limits_e {{\Omega _e}} $$ in (56) as a sum of integrals over FEs Ωe, and after that we use expression (61):  \begin{eqnarray*} {G_{lj}} &=& \frac{1}{{{\mu _0}}}\int_{\Omega }{}\left( {\nabla \times \vec{\boldsymbol\psi }_l^c} \right) \cdot \left( {\nabla \times \vec{\boldsymbol\psi }_j^c} \right){\rm{d}}\Omega = \frac{1}{{{\mu _0}}}\sum\limits_e {\int_{{{\Omega _e}}}{}\left( {\nabla \times \vec{\boldsymbol\psi }_l^c} \right) \cdot \left( {\nabla \times \vec{\boldsymbol\psi }_j^c} \right){\rm{d}}\Omega }\nonumber \\ &=& \frac{1}{{{\mu _0}}}\sum\limits_e {\int_{{{\Omega _e}}}{}\left( {\nabla \times \sum\limits_k {{T_{lk}}\vec{\boldsymbol\psi }_k^{nc}} } \right) \cdot \left( {\nabla \times \sum\limits_m {{T_{jm}}\vec{\boldsymbol\psi }_m^{nc}} } \right){\rm{d}}\Omega}. \end{eqnarray*} Now the summation over k and m can be taken off the sign of integration:   \begin{equation} {G_{lj}} = \frac{1}{{{\mu _0}}}\sum\limits_e {\sum\limits_k {\sum\limits_m {{T_{lk}}{T_{jm}}} } \int_{{{\Omega _e}}}{}\left( {\nabla \times \vec{\boldsymbol\psi }_k^{nc}} \right) \cdot \left( {\nabla \times \vec{\boldsymbol\psi }_m^{nc}} \right){\rm{d}}\Omega } . \end{equation} (63)Obviously, the integral in (63) defines the matrix $${{{\bf \hat{G}}}^e}$$ of the FE Ωe for basis functions $$\vec{\boldsymbol\psi }_k^{nc}$$ which are the standard shape vector functions on Ωe. Thus, formula (63) defines the rule of assembling matrices $${{{\bf \hat{G}}}^e}$$ when the ‘conforming’ basis functions $$\vec{\boldsymbol\psi }_j^c$$ are used on the non-conforming meshes, and shape functions on FEs (and, consequently, their matrices) are defined in a standard way for the vector FEM. This rule is as follows. The numbers k and m of global basis functions from the basis $$\{ {\vec{\boldsymbol\psi }_i^{nc}} \}$$ are defined by the numbers of local basis functions, and the considered component of the local matrix is added with the factor TlkTjm only to those components Glj of the global matrix G, for which both components Tlk and Tjm are non-zero. The generation of matrices Mσ and Mε of FE equation system (56) is performed similarly, that is, with the use of standard shape vector functions, which correspond to the ‘non-conforming’ basis $$\{ {\vec{\boldsymbol\psi }_k^{nc}} \}$$. Because the basis functions $$\vec{\boldsymbol\psi }_j^c$$ are the elements of the space Hcurl(Ω), all theoretical estimations of the convergence of FE solutions obtained with the use of these functions (Monk 2003) are true. If we perform the local mesh refinements in the needed subdomains (in the vicinity of the receivers and anomalies with significant influence), FE solutions obtained on such non-conforming meshes have essentially the same accuracy as FE solutions obtained on regular meshes, which correspond to these irregular non-conforming meshes. However, when such irregular meshes are used, the computational costs are lowered by an order or more (frequently, even more than 100 times) in comparison with the use of corresponding (i.e. with the cells refinement in the same subdomains) regular meshes. This will be shown below, in Section 5. 4.3 The grouping of logging tool positions and using the direct solver The application of the primary-secondary field approach with a separate calculation of the horizontally layered medium field and application of the unstructured optimized meshes with local refinements of the cells for calculating the secondary field allow us to significantly reduce the size of matrix of the FE system. In this case, the direct solvers become sufficiently effective. Because the solution of logging problems implies the computation of the large number of fields in the same medium for the different positions of transmitters and receivers, the application efficiency of the direct solvers can be highly increased. If the medium does not change and only the transmitter position changes, then the matrix of the FE system will remain the same for all these problems when the same mesh is used. In this case, only the right-hand side vectors of the FE system will be different (in accordance with the transmitter position). Thus, when the direct solvers for SLAE are used, the matrix factorization into triangular matrices can be performed only once for all transmitter positions, and when the field for each transmitter position is computed, one should solve only the systems with these triangular matrices (Da Silva et al.2012; Oldenburg et al.2013; Haber & Schwarzbach 2014). Note that in our code, we use Intel MKL PARDISO solver (Schenk & G̈artner 2004). However, simultaneous field calculation for all logging tool positions on the trajectory is inappropriate. First, the trajectories in the logging problems usually significantly exceed the size of the computational domain, in which the field for one transmitter position can be computed with a required accuracy. Second, in order to obtain the numerical solution with the required accuracy, it will suffice to refine the mesh only in a reasonably small vicinity of the tool position. Therefore, in order to achieve the required accuracy of the solution on the same mesh for several logging tool positions, it is required to refine the mesh in quite a large subdomain of the computational domain, and the number of nodes in the mesh (and, consequently, the SLAE size) can increase very much with the increase in the number of tool positions for which the same matrix of the FE system will be generated. As a result, the size of the FE system can increase so much that even the time of the solution of two systems with two triangular matrices (i.e. factors obtained due to the factorization) can exceed the total time of calculating the field for one tool position on the mesh with a small number of FEs (not to mention other computational costs: factorization of the big size matrix, recalculation of the primary field onto a fine mesh, and etc.). Thus, both the approach in which the field calculation is performed for each logging tool position on the individual small mesh and the approach in which the fields for all logging tool positions on the trajectory are calculated simultaneously (on the corresponding mesh with a great number of elements) are not optimal. Therefore, the most appropriate is the grouping of problems, which correspond to different tool positions on the trajectory, such that, on the one hand, there are more positions in one group, and, on the other hand, the size of the matrix of the FE system does not increase too much, because with the increase in the size of SLAE, the computational costs for matrix factorization (in order to construct the triangular factors) and other overhead costs grow very fast. In the case of increasing the number of positions in the group, the difference between the 3-D medium and the horizontally layered background medium taken for this group can also increase. In this case, in order to achieve the required accuracy, it is necessary to use a more fine mesh for calculating the secondary field that can lead to higher computational costs. We also note that the field calculation for groups can be performed in parallel that makes it possible to achieve the additional acceleration. In our code, the grouping of tool positions is performed automatically. In this case, the number of positions in a group is essentially defined by the number of positions getting into the sequent segment of the trajectory. The length of the trajectory segment (as well as the step sizes in the mesh) depends on the frequency. Thus, for the frequency of 100 kHz, the optimal length of the trajectory segment varies from 2 m (in the case, if the distance between adjacent positions is about 0.1 m) to 6 m (in the case, if the distance between adjacent positions is about 1–2 m). For other frequencies, the length of the trajectory segment is scaled by multiplying by the square root of the frequencies ratio. Thus, the number of positions in one group varies on the average from 3 to 20. 4.4 Automatic choice of the horizontally layered background medium for calculating the primary field The highest accuracy of the FE solution can be achieved without increasing the computational costs in the situation when the logging tool (i.e. its transmitter and receivers) is entirely located in such a subdomain of the 3-D computational domain the conductivity of which coincides with the conductivity of the background medium (i.e. the horizontally layered medium taken as a background for calculating the primary field). Moreover, the structures of geoelectrical models and logging tool motion trajectories, which are specific for the problems considered here (see e.g. Fig. 3), as a rule, make practically impossible the choice (pre-setting) of the single (for the entire trajectory) background medium, which would satisfy the requirements mentioned above. Therefore, we propose to perform the automatic choice of the background medium for each logging tool position (or group of positions) on the trajectory taking into account the location of the tool and the properties of the medium in close vicinity of it. For calculating the field for one group of logging tool positions, we propose the following algorithm of the automatic choice of the background medium. In the 3-D FE mesh, we select the subdomain ΩT-R, which in plane contains all transmitters and receivers of the considered group and some vicinity surrounding them (i.e. this subdomain is limited only in x- and y-directions, but in the z-direction, it coincides with the original 3-D computational domain). Using all FEs, which are inside this subdomain, the array ZT-R of z-coordinates of these FEs boundaries is generated. This array is ordered, and the numbers of FEs, the z-coordinates of which are located inside the interval ($$Z_i^{{\rm{T - R}}}$$, $$Z_{i + 1}^{{\rm{T - R}}}$$), are assigned to this interval. This information is formed as a list by one looping over all FEs from ΩT-R. Then, by means of looping over intervals the boundaries of which are defined by the array ZT-R, the horizontally layered (background) medium is formed by the following algorithm. For the current interval, we define the maximal conductivity in this interval by analysing all FEs associated with this interval and taken from the corresponding list. If this conductivity coincides with the maximal conductivity obtained for the previous interval, then we pass to the next interval. If we obtain a different value of maximal conductivity, then the interface between the current and previous intervals becomes the interface of the horizontally layered (background) medium to be formed. As conductivities of this horizontally layered (background) medium we take maximal conductivity values, which are obtained during the determination of the z-coordinates of layers interfaces. Because the size of the 3-D computational domain is smaller than the size of the computational domain for calculating the primary field, outside ΩT-R, the background horizontally layered medium is formed using the data describing the original problem. At the end of the algorithm, we check the thickness of the layers. If the layers with a very small thickness are presented in the formed medium, then these layers are joined with the adjacent ones (from two adjacent layers we choose the layer which is thicker), that is, in fact, a small thickness layer is excluded from the medium, and the layer adjacent to it is properly extended. We note that the 3-D mesh, which is used for the automatic choice of the background medium, can partially describe the geometry of the 3-D problem to be solved. Before automatically choosing, it is better not to include small size objects into this mesh (e.g. a borehole if its effect on the measured signal is rather small). We include these objects into the 3-D mesh after the layers of the background horizontally layered medium have been formed, and their interfaces have been included into the mesh. Fig. 9 (right) shows the sections of the background mediums, which are formed with the use of the algorithm considered above for three groups of tool positions, including the situation when the trajectory crosses the fault. Figure 9. View largeDownload slide The vertical sections of meshes (left) generated for different groups of adjacent tool positions for the case when the background medium (right) is chosen automatically. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the ith transmitter. Figure 9. View largeDownload slide The vertical sections of meshes (left) generated for different groups of adjacent tool positions for the case when the background medium (right) is chosen automatically. Ti denotes the location of the ith transmitter; Ri,j denotes the location of the jth receiver of the ith transmitter. Fig. 9 shows the FE meshes generated for the same tool positions on the trajectory as shown in Figs 5 and 6. But in Fig. 9, these meshes have been generated for the case when the background medium is chosen automatically for each group of the tool positions. The comparison of the numerical results for the cases when the background medium is pre-set as single and when the background medium is chosen automatically (for each group of the adjacent tool positions) will be presented below, in Section 6. 5 VALIDATION In this section, we present the results of the verification of the mathematical models and numerical schemes proposed in this paper and described in the previous sections. We will call the field calculation performed with the use of mathematical models (3), (20), (28)–(30) the 2-D approximation because, although the medium is horizontally layered (i.e. really 1-D), in order to calculate the electromagnetic field induced by the current coil, we use the solution of one or several axisymmetrical (i.e. 2-D) problems in the 2-D computational domain in the (r, z)-coordinates. We will call the field calculation performed with the use of mathematical model (1) the 3-D approximation even in the case when the 2.5-D problem is solved, because in order to obtain its solution, we find the FE solution in the 3-D computational domain. In order to verify the method proposed in this paper for calculating the primary field in the horizontally layered medium, we use the Dipole1D code (Key 2009). Because the Dipole1D code allows the calculation of the field induced only by the electrical dipole, then in order to obtain the field induced by the current frame (loop), we sum the fields of four dipoles making up the current frame. As a testing model, we use the geoelectrical model shown in Fig. 10. Figure 10. View largeDownload slide Testing geoelectrical model 1 of the horizontally layered medium and locations of the lines (Line1…Line5) along which the non-zero components of $${{\bf \vec{B}}}$$ and $${{\bf \vec{E}}}$$ are calculated with the use of the proposed approach and Dipole1D code (Key 2009). Figure 10. View largeDownload slide Testing geoelectrical model 1 of the horizontally layered medium and locations of the lines (Line1…Line5) along which the non-zero components of $${{\bf \vec{B}}}$$ and $${{\bf \vec{E}}}$$ are calculated with the use of the proposed approach and Dipole1D code (Key 2009). The source in the form of a square frame with the side of 0.1 m is fixed, and its centre is located at the point (0, 0, 13). The current is 100 A. The frequency is 100 kHz. In order to verify the field values in the receivers, we will output the non-zero components of the magnetic induction $${{{\bf \vec{B}}}^p}$$ along Line 1 and Line 2 (see Fig. 10). Besides, as we have described in the previous sections, when calculating the 3-D field, the right-hand side of the equations (see eq. 1) is defined by the electrical field induced by the source in the horizontally layered medium, which has been taken as a background. Therefore, we will also verify the accuracy of calculating the non-zero components of the electrical field $${{{\bf \vec{E}}}^p}$$. We will analyse the accuracy of calculating $${{{\bf \vec{E}}}^p}$$ along Line 3, Line 4, and Line 5 (see Fig. 10), and also along Line 6 and Line 7, which are parallel to Line 3 and have the same z-coordinate. The distance along the y-axis between these lines (Line 6 and Line 7) and Line 3 is 2 m and 4 m, respectively. The plots placed at the top of Fig. 11 show the real and imaginary x- and z-components of $${{{\bf \vec{B}}}^p}$$ along Line 1 and Line 2 for the case when the field is induced by the ‘vertical’ frame and ‘horizontal’ loop. The plots placed at the bottom of Fig. 11 show the relative (in per cent) deviations of $${{{\bf \vec{B}}}^p}$$, which are calculated using the approach proposed in this paper, from the values of $${{{\bf \vec{B}}}^p}$$, which are calculated with the use of the Dipole1D code (Key 2009). The results presented in Fig. 11 show that within the distance of 10 m, the relative deviations for all components do not exceed 0.5 per cent, and within the distance of 20 m, the deviations do not exceed 1–1.5 per cent. The error growth takes place in the subintervals of low values of the field. We do not consider the field values farther than 20 m, because this value (20 m) defines the maximal size of the 3-D computational domain. Figure 11. View largeDownload slide $$B_x^{{\rm{Re}}}$$ (dark blue curve), $$B_x^{{\mathop{\rm Im}\nolimits} }$$ (green curve), $$B_z^{{\rm{Re}}}$$ (red curve) and $$B_z^{{\mathop{\rm Im}\noli