CSEM forward modelling & Edge elements formulation

The last decade has been a period of rapid growth for electromagnetic methods (EM) in geophysics, mostly because of their industrial adoption. In particular, the 3D marine controlled-source electromagnetic (3D CSEM) method has become an important technique for reducing ambiguities in data interpretation in hydrocarbon exploration. In order to be able to predict the EM signature of a given geological structure, modelling tools provide us with synthetic results which we can then compare to measured data. In particular, if the geology is structurally complex, one might need to use methods able to cope with such complexity in a natural way by means of, e.g., an unstructured mesh representing its geometry. Among the modelling methods for EM based upon 3D unstructured meshes, the Nédélec Edge Finite Element Method (EFEM) offers a good trade-off between accuracy and number of degrees of freedom, e.g. size of the problem. Furthermore, its divergence-free basis is very well suited for solving Maxwell’s equation. On top of that, we choose to support tetrahedral meshes, as these are the easiest to use for very large domains or complex geometries. We present the numerical formulation and results of 3D CSEM forward modelling (FM) using EFEM of high-order on unstructured tetrahedral meshes.

CSEM problem

3D CSEM FM is typically solved in frequency domain, which involves the numerical solution of Maxwell’s equations in stationary regimes for heterogeneous anisotropic electrically conductive domains. Furthermore, CSEM surveys generally work with low frequency electromagnetic fields (\sim 1 Hz to \sim 3 Hz) because the electric conductivity of the geological structures is much larger than their dielectric permittivity. As a consequence, in an unbound domain \Gamma, the electric field can be obtained by solving Maxwell’s equations in their diffusive form

(1)\nabla \times \mathbf{E} &= i\omega \mu_{0}\mathbf H,  \\
      \nabla \times \mathbf{H} &= \mathbf J_{s} + \tilde{\sigma}\mathbf E,

where the harmonic time dependence e^{-i \omega t} is omitted, with \omega is the angular frequency, \mu_{0} the free space magnetic permeability, \mathbf J_{s} the distribution of source current, \tilde{\sigma}\mathbf E the induced current in the conductive Earth and \tilde{\sigma} the electrical conductivity which is assumed isotropic for simplicity.

In numerical approximations of EM fields there are two main drawbacks. The first one is the inevitable spatial singularity at the source. The second is the grid refinement requirements in order to capture the rapid change of the primary field [1]. In order to mitigate these issues, PETGEM used a secondary field approach where the total electric field \mathbf E is obtained as

(2)\mathbf E &= \mathbf E_{p} + \mathbf E_{s}, \\
\tilde{\sigma} &=  \tilde{\sigma_{s}} + \Delta \tilde{\sigma},

where subscripts p and s represent a primary field and secondary field respectively. For a general layered Earth model, \mathbf E_{p} can be computed semi-analytically by using Hankel transform filters. Based on this decomposition and following the work by [3] the equation system to solve \mathbf E_{s} is

(3)\nabla \times \nabla \times \mathbf E_{s} + i \omega \mu \tilde{\sigma} \mathbf E_{s} = -i \omega \mu \Delta \sigma \mathbf E_{p},

where the electrical conductivity \sigma is a function of position that is allowed to vary in 3D, whereas the vacuum permeability \mu is set to the free space value \mu_{0}. Homogeneous Dirichlet boundary conditions, \mathbf E_{s} = 0 on \partial\Gamma, are applied. The range of applicability of this conditions can be determined based on the skin depth of the electric field [4].

Edge finite element formulation

For the computation of \mathbf E_{s}, PETGEM implemented the high-order Nédélec EFEM which uses vector basis functions defined on the edges, faces and volume of the corresponding elements. Its basis functions are divergence-free but not curl-free [2]. Thus, EFEM naturally ensures tangential continuity and allows normal discontinuity of \mathbf E_{s} at material interfaces. PETGEM used unstructured tetrahedral meshes because of their ability to represent complex geological structures such as bathymetry or reservoirs as well as the local refinement capability in order to improve the solution accuracy. Figure 4.1 shows the tetrahedral Nédélec elements (first-, second-, and third-order) together with their node and edge indexing.

Tetrahedral Nédélec edge element with node/edge/face indexing.

Figure 4.1. Reference tetrahedral element showing numeration of DOFs for each order p. For p=1 there are 6 DOFs (1 per edge). For p=2 there are 20 DOFs (2 per edge and 2 per face). Finally, for p=3 there are 45 DOFs (3 per edge, 6 per face and 3 per element’s volume).

In PETGEM workflow, the tangential component of the secondary electric field is assigned to the edges in the mesh. Therefore, all components of the electric field at a point \mathbf x located inside a tetrahedral element e can be obtained as follows

(4)\mathbf E^{e}(\mathbf x) = \sum_{i=1}^{ndofs} \mathbf N^{e}_{i}(\mathbf x) E^{e}_{i},

where \mathbf N^{e}_{i} are the vector basis functions and E^{e}_{i} their degrees of freedom. For instance, considering the node and edge indexing in Figure 4.1, the vector basis functions for p=1 (first-order) can be expressed as follows

(5)\mathbf N^{e}_{i} &= (\lambda^{e}_{i1} \nabla \lambda^{e}_{i2} - \lambda^{e}_{i2} \nabla \lambda^{e}_{i1}) \ell^{e}_{i},

where subscripts i1 and i2 are the first and second nodes linked to the i-th edge, \lambda^{e}_{i} are the linear nodal basis functions, and \ell^{e}_{i} is the length of the i-th edge of the element e. A systematic approach for obtaining mixed-order curl-conforming basis functions may be seen in [4], [5].

By substituting equation (4) into (3), and using Galerkin’s approach, the weak form of the original differential equation becomes

(6)Q_{i} = \int_{\Omega} \mathbf N_{i} \cdot [ \nabla \times \nabla \times \mathbf E_{s} -i \omega \mu \tilde{\sigma} \mathbf E_{s} + i \omega \mu \Delta \tilde{\sigma} \mathbf E_{p} ] dV,

The compact discretized form of (6) is obtained after applying the Green’s theorem

(7)[K^{e}_{jk} + i \omega \tilde{\sigma}_{e} M^{e}_{jk}] \cdot  \{ E_{sk} \} = - i \omega \mu \Delta \tilde{\sigma}_{e} R^{e}_k,

where K^{e} and M^{e} are the elemental stiffness and mass matrices which can be calculated analytically or numerically [2], and R^{e}_k is the right hand side which requires numerical integration.

[1]Cai, H., Xiong, B., Han, M. and Zhdanov, M. (2014). 3D controlled-source electromagnetic modeling in anisotropic medium using edge-based finite element method. Computers & Geosciences, 73, 164–176.
[2](1, 2) Jin, J. (2002). The Finite Element Method in Electromagnetics. Wiley, New York, second edn.
[3]Newman, G.A. and Alumbaugh, D.L. (2002). Three-dimensional induction logging problems, Part 2: A finite difference solution. Geophysics, 67(2), 484–491.
[4](1, 2) Puzyrev, V., Koldan, J., de la Puente, J., Houzeaux, G., Vázquez, M. and Cela, J.M. (2013). A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling. Geophysical Journal International, ggt027.
[5]Garcia-Castillo, L.E., Salazar-Palma, M., 2000. Second-order Nédélec tetrahedral element for computational electromagnetics. International Journal of Numerical Modelling: Electronic Networks, Devices and Fields (John Wiley & Sons, Inc.) 13, 261–287.
[6]Garcia-Castillo, L.E., Ruiz-Genovés, A.J., Gómez-Revuelto, I., Salazar-Palma, M., Sarkar, T.K., 2002. Third-order Nédélec curl-conforming finite element. IEEE Transactions on Magnetics 38, 2370–2372.