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 marine controlled-source electromagnetic (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 real 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, i.e. 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 tetrahedral EFEM on unstructured 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) 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)\[\begin{split}\nabla \times \mathbf{E} &= i\omega \mu_{0}\mathbf H \\ \nabla \times \mathbf{H} &= \mathbf J_{s} + \tilde{\sigma}\mathbf E\end{split}\]

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)\[\begin{split}\mathbf E &= \mathbf E_{p} + \mathbf E_{s} \\ \tilde{\sigma} &= \tilde{\sigma_{s}} + \Delta \tilde{\sigma}\end{split}\]

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 [6] 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 [7].

Edge finite element formulation

For the computation of \(\mathbf E_{s}\), PETGEM implemented the Nédélec EFEM which uses vector basis functions defined on the edges of the corresponding elements. Its basis functions are divergence-free but not curl-free [4]. 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 (lowest order) together with their node and edge indexing.

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

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

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}^{6} \mathbf N^{e}_{i}(\mathbf x) E^{e}_{i}\]

where \(\mathbf N^{e}_{i}\) are the vector basis functions associated to each edge \(i\) and \(E^{e}_{i}\) their degrees of freedom. Considering the node and edge indexing in Figure 4.1, the vector basis functions can be expressed as follows:

(5)\[\begin{split}\mathbf N^{e}_{i} &= (\lambda^{e}_{i1} \nabla \lambda^{e}_{i2} - \lambda^{e}_{i2} \nabla \lambda^{e}_{i1}) \ell^{e}_{i}\end{split}\]

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\).

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 [4], 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]Castillo, O., de la Puente, J., Puzyrev, V. and Cela, J.M. (2015). Edge-based electric field formulation in 3D CSEM simulations: a parallel approach. In: Proceedings of the 6th International Conference and Workshop on Computing and Communication. IEEE. Vancouver, Canada.
[3]Constable, S. and Weiss, C.J. (2006). Mapping thin resistors and hydrocarbons with marine EM methods:Insights from 1D modeling. Geophysics, 71(2), G43–G51.
[4](1, 2) Jin, J. (2002). The Finite Element Method in Electromagnetics. Wiley, New York, second edn.
[5]Key, K. (2009). 1D inversion of multicomponent, multifrequency marine CSEM data: Methodology and synthetic studies for resolving thin resistive layers, F9–F20.
[6]Newman, G.A. and Alumbaugh, D.L. (2002). Three-dimensional induction logging problems, Part 2: A finite difference solution. Geophysics, 67(2), 484–491.
[7]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.