A New 3D Marine Controlled-Source Electromagnetic Modeling Algorithm Based on Two New Types of Quadratic Edge Elements

Abstract

This study proposes a new algorithm for a higher-order vector finite element method based on two new types of second-order edge elements to solve the electromagnetic field diffusion problem in a 3D anisotropic medium. To avoid source singularity in the quasistatic variant of Maxwell’s function, a secondary field formulation was adopted. The modeling domain was discretized using two types of quadratic edge hexahedral elements, which were obtained using the edge unification method to reduce variables on each side of two conventional quadratic edge elements. Compared with the traditional quadratic element, the number of unknowns that needed to be solved was significantly reduced. The sparse linear equation of the finite element system was solved using an open-source direct solver called MUMPS. The numerical results demonstrated that the proposed algorithm has the same level of accuracy as the conventional vector finite element method and has a significant advantage over it in terms of computational cost.

1. Introduction

The marine controlled-source electromagnetic method (MCSEM) uses controllable artificial field sources to transmit electromagnetic signals in the ocean and measures the electric or magnetic fields at a location far away from the field source to detect the electrical distribution of the media below the seabed. The resistance characteristics of oil and gas reservoirs are the most important properties, which can generate characteristic surface electromagnetic signals. In other words, MCSEM technology can distinguish underground oil and gas from other fluids. In theory, CSEM measurement is to record data with multiple power-receiver offsets, several different frequencies, and a certain power-receiver arrangement. The marine controlled-source electromagnetic (MCSEM) method has become a popular geophysical exploration tool for offshore hydrocarbon (HC) exploration [1,2,3]. During the last decade, marine CSEM has been widely used to reduce ambiguities in data interpretation and reduce the risk of exploration. With the rapid development of MCSEM exploration instruments, large-scale and high-quality data have been obtained, resulting in the need to develop high-precision MCSEM three-dimensional (3D) data interpretation tools. The 3D inversion algorithm technology of MCSEM is commonly used for data interpretation. Forward modeling is the basis of inversion, and high-precision forward modeling technology is key to achieving high-precision inversion.

The integral equation (IE) [4,5,6], finite difference (FD) [7,8,9,10,11,12], and finite element (FE) methods [13,14,15,16,17] are popular numerical techniques for electromagnetic (EM) field modeling. In the FE method, it is easier to use an unstructured mesh to create an irregular body, locally refine the interest region, and coarsen the boundary area of the model domain [15]. With this processing, numerical modeling becomes efficient and effective. Therefore, these advantages make the FE method more suitable for simulating the complex EM response of irregular models than FD and IE methods.

Compared with the conventional finite element method, which is based on a node, the edge-based finite element method has the advantage that divergence-free conditions are automatically satisfied by appropriately selected basic functions [18]. In recent years, the vector finite element method based on edge elements has been used to solve the electromagnetic field problem, and is known as Maxwell’s equations [19,20]. It should be noted that to date, most formulations for the 3D CSEM problem have implemented the lowest order. In addition, some studies have used higher-order elements to solve three-dimensional forward modeling of CSEM, but these were all aimed at tetrahedral elements [14,21]. To the best of our knowledge, no high-order hexahedron element is currently used for three-dimensional forward modeling of the controlled-source electromagnetic method. As we know, high-order edge elements can provide accurate solutions but result in an increase in unknown variables and lead to larger computing costs. Fortunately, we can improve the accuracy of the solution at a small computational cost by eliminating some variables when using higher-order elements [22,23].

In this study, we applied the unification method to the edge to reduce one edge variable on each side of the conventional serendipity and the Lagrange-type quadratic hexahedral element, resulting in two new types of second order edge elements. We also introduced a new vector FE method based on the proposed elements to solve the 3D marine CSEM modeling problem in an anisotropic medium. For a complex bathymetry simulation, we used a general hexahedral element to discretize the modeling domain. We validated our code by using several models in isotropic and anisotropic media.

2. Theory

In geophysical applications, the displacement current is typically ignored when the frequency of the electromagnetic field is low. Maxwell’s equations can then be simplified as follows [24]:

$$\nabla \times \mathit{E}=i\omega {\mu }_0\mathit{H}$$

(1)

$$\nabla \times \mathit{H}=\sigma \mathit{E}+{\mathit{J}}_{\mathrm{s}}$$

(2)

where the space magnetic permeability is denoted by ${\mu }_0$, angular frequency is denoted by $\omega$, source current is denoted by ${\mathit{J}}_{\mathrm{s}}$, and conductivity tensor is denoted by $\sigma$:

$$\mathit{\sigma }=\left(\begin{array}{ccc}{\sigma }_x& 0& 0\\ 0& {\sigma }_y& 0\\ 0& 0& {\sigma }_z\end{array}\right)$$

(3)

where ${\sigma }_x$, ${\sigma }_y$ and ${\sigma }_z$ are principal conductivities. In this study, we discuss only the principal axis anisotropy case. However, our formulation also works in the case of general anisotropy.

Calculate the curl at the left and right ends of Equation (1), and use Equation (2) to eliminate the magnetic field:

$$\nabla \times \nabla \times \mathit{E}-i\omega {\mu }_0\mathit{\sigma }\mathit{E}=i\omega {\mu }_0{\mathit{J}}_{\mathrm{s}}$$

(4)

The total and secondary field formulations are two types of formulations that are typically used for solving the 3D MCSEM modeling problem. In contrast to the total field formulation, the secondary field formulation’s source term comprises a primary electric field. This type of source term is smoother than the current source used in the total-field formulation. Therefore, in this study, we use a secondary field formulation to address the marine CSEM modeling problem in an anisotropic medium. In this formulation, the total field can be decomposed into two parts: the background (primary field, ${\mathit{E}}_{\mathrm{p}}$) and the anomaly fields (secondary field, ${\mathit{E}}_{\mathrm{s}}$):

$$\mathit{E}={\mathit{E}}_{\mathrm{s}}+{\mathit{E}}_{\mathrm{p}}$$

(5)

Similarly, the conductivity tensor can be decomposed into background (${\mathit{\sigma }}_{\mathrm{b}}$) and anomalous conductivities ($\mathsf{\Delta }\mathit{\sigma }$):

$$\mathit{\sigma }={\mathit{\sigma }}_{\mathrm{b}}+\mathsf{\Delta }\mathit{\sigma }$$

(6)

With substitution of Equations (5) and (6) into Equation (4), the Helmholtz equation with the electric component of the secondary field can be written as:

$$\nabla \times \nabla \times {\mathit{E}}_{\mathrm{s}}-i\omega {\mu }_0\mathit{\sigma }{\mathit{E}}_{\mathrm{s}}=i\omega {\mu }_0\mathsf{\Delta }\mathit{\sigma }{\mathit{E}}_{\mathrm{p}}$$

(7)

The secondary magnetic component can easily be calculated using Equation (8) once the secondary electric component in Equation (7) can be solved:

$${\mathit{H}}_{\mathrm{s}}=\frac{1}{i\omega {\mu }_0}\nabla \times {\mathit{E}}_{\mathrm{s}}$$

(8)

3. Edge-Based Finite Element Analysis

In this study, we used the edge-based FE method with quadratic elements, as discussed in the previous section. We used the edge shape function ${\mathit{N}}_i$ to approximate the secondary electric field within an edge element:

$${\mathit{E}}_s={\displaystyle \sum }_i^{N_{\mathrm{edge}}}{\mathit{N}}_i{\mathit{E}}_{s,j}$$

(9)

where $N_{\mathrm{edge}}$ denotes the total number of edges in the edge element. For the first-order edge element, the $N_{\mathrm{edge}}=12$; for conventional serendipity and Lagrange-type second-order hexahedral elements, $N_{\mathrm{edge}}=36$ and $54$, respectively.

By substituting (9) into (7) and applying the Galerkin method, we can determine that the weak form of the original differential equation is:

$${\mathrm{R}}_i={{\displaystyle \int }}_{\mathsf{\Omega }}{\mathit{N}}_i·\left[\nabla \times \nabla \times {\mathit{E}}_s-i\omega {\mu }_0\mathit{\sigma }{\mathit{E}}_s-i\omega {\mu }_0\mathsf{\Delta }\mathit{\sigma }{\mathit{E}}_{\mathrm{p}}\right]\mathrm{d}v$$

(10)

where $\mathsf{\Omega }$ denotes the modeling area.

Applying first vector Green’s theorem to Equation (10), the secondary electric field is given, and it is continuous across the inner boundary; the surface term in Equation (10) vanishes. Then, the discretized form of Equation (10) for each element can be expressed as follows:

$${\mathrm{R}}_i^{\mathrm{N}e}={\displaystyle \sum }_1^{\mathrm{N}e}\left[{\mathit{K}}^e{\mathit{E}}_{\mathrm{s}}^{\mathrm{e}}-i\omega {\mu }_0{\mathit{M}}_1^e{\mathit{E}}_{\mathrm{s}}^e-i\omega {\mu }_0{\mathit{M}}_2^e{\mathit{E}}_{\mathrm{p}}^e\right]$$

(11)

where ${\mathit{K}}^e$ and ${\mathit{M}}_1^e$ and ${\mathit{M}}_2^e$ are the local stiffness matrices defined as follows:

$${\mathit{K}}_{i,j}^e={{\displaystyle \int }}_{{\mathsf{\Omega }}_e}\left(\nabla \times {\mathit{N}}_i^e\right)·\left(\nabla \times {\mathit{N}}_j^e\right)\mathrm{d}v$$

(12)

$${\mathit{M}}_{1i,j}^{\mathrm{e}}={{\displaystyle \int }}_{{\mathsf{\Omega }}_{\mathrm{e}}}{\mathit{N}}_i^{\mathrm{e}}·\mathit{\sigma }·{\mathit{N}}_j^{\mathrm{e}}\mathrm{d}v$$

(13)

$${\mathit{M}}_{2i,j}^{\mathrm{e}}={{\displaystyle \int }}_{{\mathsf{\Omega }}_{\mathrm{e}}}{\mathit{N}}_i^{\mathrm{e}}·\mathsf{\Delta }\mathit{\sigma }·{\mathit{N}}_j^{\mathrm{e}}\mathrm{d}v$$

(14)

where ${\mathsf{\Omega }}_{\mathrm{e}}$ denotes the domain of an element.

By assembling a global matrix with local matrices, a linear equation for the elements is obtained. The equation is typically sparse:

$$\mathit{A}\mathit{e}=\mathit{b},$$

(15)

Here $\mathit{A}$ denotes the global stiffness matrix, and $\mathit{b}$ denotes the RHSs of the equations.

To solve Equation (14), a proper condition should be added to the equation. For simplicity, let us consider the homogeneous Dirichlet boundary condition:

$$\mathit{e}{|}_{\mathsf{\Gamma }}=0$$

(16)

Here, $\mathsf{\Gamma }$ is the boundary area.

IDR(s) iterative solvers with an ILU preconditioner were used to solve the linear equations in Equation (15).

4. First and Conventional Quadratic Edge Elements

In the edge-based FE method, we can use regular rectangular, general tetrahedral, general hexahedral, or other complex elements to discretize the modeling domain. To simplify the problem and simulate bathymetry, hexahedral elements were used to discretize the modeling domain in this study.

To calculate the stiffness matrix of the hexahedral element, a general hexahedral element should be transformed into an origin-centered cubic element. A general hexahedron in Cartesian coordinates $xyz$ is shown in Figure 1a, and a transformed reference domain (cubic element) in Cartesian coordinates $\xi \zeta \eta$ is shown in Figure 1b.

The transformation can be described by the following formulas [25]:

$$x={{\displaystyle \sum }}_{i=1}^8{N}_i^e\left(\xi ,\eta ,\zeta \right)x_i^e$$

(17)

$$y={{\displaystyle \sum }}_{i=1}^8{N}_i^e\left(\xi ,\eta ,\zeta \right)y_i^e$$

(18)

$$z={{\displaystyle \sum }}_{i=1}^8{N}_i^e\left(\xi ,\eta ,\zeta \right)z_i^e$$

(19)

where $N_i^e$ is the scalar node-based shape function, which can be defined as follows:

$$N_i^e\left(\xi ,\eta ,\zeta \right)=\frac{1}{8}\left(1+\xi {\xi }_i\right)\left(1+\eta {\eta }_i\right)\left(1+\zeta {\zeta }_i\right)$$

(20)

where $i$ is a local node index of the element.

In the analysis of the vector finite element, two volume integrals (12)–(14) must be evaluated. For a hexahedron element, it is easy to transform with two integrals in the $\xi \eta \zeta$-coordinate using the Jacobian transform:

$$K_{i,j}^e={{\displaystyle \int }}_{-1}^1{{\displaystyle \int }}_{-1}^1{{\displaystyle \int }}_{-1}^1\left(\nabla \times {\mathit{N}}_i^e\left(\xi ,\eta ,\zeta \right)\right)·\left(\nabla \times {\mathit{N}}_j^e\left(\xi ,\eta ,\zeta \right)\right)\left|\mathit{J}\left(\xi ,\eta ,\zeta \right)\right|\mathrm{d}\xi \mathrm{d}\eta \mathrm{d}\zeta ,$$

(21)

$$M_{1i,j}^e={{\displaystyle \int }}_{-1}^1{{\displaystyle \int }}_{-1}^1{{\displaystyle \int }}_{-1}^1{\mathit{N}}_i^e\left(\xi ,\eta ,\zeta \right)·\mathit{\sigma }·{\mathit{N}}_j^e\left(\xi ,\eta ,\zeta \right)\left|\mathit{J}\left(\xi ,\eta ,\zeta \right)\right|\mathrm{d}\xi \mathrm{d}\eta \mathrm{d}\zeta ,$$

(22)

$$M_{2i,j}^e={{\displaystyle \int }}_{-1}^1{{\displaystyle \int }}_{-1}^1{{\displaystyle \int }}_{-1}^1{\mathit{N}}_i^e\left(\xi ,\eta ,\zeta \right)·\mathsf{\Delta }\mathit{\sigma }·{\mathit{N}}_j^e\left(\xi ,\eta ,\zeta \right)\left|\mathit{J}\left(\xi ,\eta ,\zeta \right)\right|\mathrm{d}\xi \mathrm{d}\eta \mathrm{d}\zeta ,$$

(23)

$\mathit{J}\left(\xi ,\eta ,\zeta \right)$ and $\left|\mathit{J}\left(\xi ,\eta ,\zeta \right)\right|$ denote the Jacobian matrix and determinant of $\mathit{J}\left(\xi ,\eta ,\zeta \right)$, respectively.

$$\mathit{J}\left(\xi ,\eta ,\zeta \right)=\left[\begin{array}{ccc}\frac{\partial x}{\partial \xi }& \frac{\partial y}{\partial \xi }& \frac{\partial z}{\partial \xi }\\ \frac{\partial x}{\partial \eta }& \frac{\partial y}{\partial \eta }& \frac{\partial z}{\partial \eta }\\ \frac{\partial x}{\partial \zeta }& \frac{\partial y}{\partial \zeta }& \frac{\partial z}{\partial \zeta }\end{array}\right]$$

(24)

In the next section, we introduce several conventional first-order and quadratic hexahedral elements and two quadratic edge hexahedral elements of the new type.

4.1. First-Order Edge Hexahedral Element

A conventional first-order hexahedral element is shown in Figure 2, which has eight nodes and twelve edges. The shape function of the edge on the sides along the $\xi$-axis direction are accordingly defined as [25]:

$${\mathit{N}}_i=\frac{1}{8}\left(1+{\eta }_i\eta \right)\left(1+{\zeta }_i\zeta \right)\mathsf{\Delta }\xi$$

(25)

Similarly, the shape function of the edge on the sides along the $\eta$-axis or $\zeta$-axis direction is defined as:

$${\mathit{N}}_j=\frac{1}{8}\left(1+{\xi }_j\xi \right)\left(1+{\zeta }_j\zeta \right)\mathsf{\Delta }\eta$$

(26)

and

$${\mathit{N}}_k=\frac{1}{8}\left(1+{\xi }_k\xi \right)\left(1+{\zeta }_k\zeta \right)\mathsf{\Delta }\zeta$$

(27)

where $\xi$, $\eta$ and $\zeta$ denote the local coordinate in the elements.

4.2. A Quadratic Hexahedron Element of a Conventional Serendipity Type

Figure 3 shows a conventional serendipity quadratic hexahedron element with 12 edges on its surface, 20 nodes in the element, and 24 edges on the sides. The shape function for the serendipitous quadratic hexahedron element can be defined as follows: the shape function of the edge on sides which along $\xi$-axis direction should be:

$${\mathit{N}}_i=\frac{1}{8}\left(1+{\eta }_i\eta \right)\left(1+{\zeta }_i\zeta \right)\left(4{\xi }_i\xi +{\eta }_i\eta +{\zeta }_i\zeta -1\right)\mathsf{\Delta }\xi$$

(28)

The edge shape function of the edges along the $\xi$-axis direction on the $\zeta =\pm 1$ surface:

$${\mathit{N}}_i=\frac{1}{4}\left(1-{\eta }^2\right)\left(1+{\zeta }_i\zeta \right)\mathsf{\Delta }\xi$$

(29)

The edge shape function of the edges along the $\xi$-axis direction on the $\eta =\pm 1$ surface:

$${\mathit{N}}_i=\frac{1}{4}\left(1-{\zeta }^2\right)\left(1+{\eta }_i\eta \right)\mathsf{\Delta }\xi$$

(30)

The sides along the $\eta$-axis direction on the sides:

$${\mathit{N}}_j=\frac{1}{8}\left(1+{\xi }_j\xi \right)\left(1+{\zeta }_j\zeta \right)\left(4{\eta }_j\eta +{\xi }_j\xi +{\zeta }_j\zeta -1\right)\mathsf{\Delta }\eta$$

(31)

The edges along the $\eta$-axis direction on the $\xi =\pm 1$ surface:

$${\mathit{N}}_j=\frac{1}{4}\left(1-{\zeta }^2\right)\left(1+{\xi }_j\xi \right)\mathsf{\Delta }\eta$$

(32)

The edges along the $\eta$-axis direction on the $\zeta =\pm 1$ surface:

$${\mathit{N}}_j=\frac{1}{4}\left(1-{\xi }^2\right)\left(1+{\zeta }_j\zeta \right)\mathsf{\Delta }\eta$$

(33)

The edges along the $\zeta$-axis direction on the $\eta =\pm 1$surface:

$${\mathit{N}}_k=\frac{1}{4}\left(1-{\xi }^2\right)\left(1+{\eta }_k\eta \right)\mathsf{\Delta }\zeta$$

(34)

4.3. Lagrange Quadratic Hexahedron Element

Figure 4 shows a Lagrange quadratic hexahedron element (LC) with 27 nodes and 24 edges on the sides, 24 edges on the surface, and 8 edges within the element. The shape function of the Lagrange quadratic hexahedron element is defined as the edges on the sides along the $\xi$-axis direction:

$${\mathit{N}}_i=\frac{1}{8}\left(1+4{\xi }_i\xi \right)\eta \zeta \left({\eta }_i+\eta \right)\left({\zeta }_i+\zeta \right)\mathsf{\Delta }\xi$$

(35)

The edges along $\mathrm{the} \xi$-axis direction on the$\zeta =\pm 1$ surface:

$${\mathit{N}}_i=\frac{1}{4}\left(1+4{\xi }_i\xi \right)\zeta \left(1-{\eta }^2\right)\left({\zeta }_i+\zeta \right)\mathsf{\Delta }\xi$$

(36)

The edges along the $\xi$-axis direction on the $\eta =\pm 1$ surface:

$${\mathit{N}}_i=\frac{1}{4}\left(1+4{\xi }_i\xi \right)\eta \left(1-{\zeta }^2\right)\left({\eta }_i+\eta \right)\mathsf{\Delta }\xi$$

(37)

The edges along $\mathrm{the} \xi$-axis direction within the element:

$${\mathit{N}}_i=\frac{1}{2}\left(1+4{\xi }_i\xi \right)\left(1-{\eta }^2\right)\left(1-{\zeta }^2\right)\mathsf{\Delta }\xi$$

(38)

The edges along $\mathrm{the} \eta$-axis direction on the sides:

$${\mathit{N}}_j=\frac{1}{8}\left(1+4{\eta }_j\eta \right)\xi \zeta \left({\xi }_j+\xi \right)\left({\zeta }_j+\zeta \right)\mathsf{\Delta }\eta$$

(39)

The edges along the $\eta$-axis direction on the $\zeta =\pm 1$ surface:

$${\mathit{N}}_j=\frac{1}{4}\left(1+4{\eta }_j\eta \right)\zeta \left(1-{\xi }^2\right)\left({\zeta }_j+\zeta \right)\mathsf{\Delta }\eta$$

(40)

The edges along the $\eta$-axis direction on the $\xi =\pm 1$ surface:

$${\mathit{N}}_j=\frac{1}{4}\left(1+4{\eta }_j\eta \right)\xi \left(1-{\zeta }^2\right)\left({\xi }_j+\xi \right)\mathsf{\Delta }\eta$$

(41)

The edges along the $\eta$-axis direction within the element:

$${\mathit{N}}_j=\frac{1}{2}\left(1+4{\eta }_j\eta \right)\left(1-{\xi }^2\right)\left(1-{\zeta }^2\right)\mathsf{\Delta }\eta$$

(42)

The edges along the $\zeta$-axis direction on the sides:

$${\mathit{N}}_k=\frac{1}{8}\left(1+4{\zeta }_k\zeta \right)\xi \eta \left({\xi }_k+\xi \right)\left({\eta }_k+\eta \right)\mathsf{\Delta }\zeta$$

(43)

The edges along the $\zeta$-axis direction on the $\xi =\pm 1$ surface:

$${\mathit{N}}_k=\frac{1}{4}\left(1+4{\zeta }_k\zeta \right)\xi \left(1-{\eta }^2\right)\left({\xi }_k+\xi \right)\mathsf{\Delta }\zeta$$

(44)

The edges along the $\zeta$-axis direction on the $\eta =\pm 1$ surface:

$${\mathit{N}}_k=\frac{1}{4}\left(1+4{\zeta }_k\zeta \right)\eta \left(1-{\xi }^2\right)\left({\eta }_k+\eta \right)\mathsf{\Delta }\zeta$$

(45)

The edges along the $\zeta$-axis direction within the element:

$${\mathit{N}}_k=\frac{1}{2}\left(1+4{\zeta }_k\zeta \right)\left(1-{\eta }^2\right)\left(1-{\xi }^2\right)\mathsf{\Delta }\zeta$$

(46)

4.4. A New Type of Quadratic Edge Element

There are two edges on the side of a conventional quadratic hexahedron element. Within an edge element, because the curls of the edge shape functions are linearly dependent, many edge variables are redundant. The computation time and storage requirements of these elements has increased significantly. Therefore, eliminating these redundant variables helps improve efficiency at no cost to computational accuracy. In this section, we adopt a new method to eliminate redundant variables from the conventional quadratic element.

Let us first consider a one-sided conventional quadratic edge element, as shown in Figure 3 and Figure 4. We named the two edges on this side edge1 and edge2. ${\mathit{N}}_1$ and ${\mathit{N}}_2$ are the shape functions of edge1 and edge2, and their corresponding edge variables are $A_1$ and $A_2$. Because one of these is redundant, we decided to eliminate $A_2$ and consider $A_1=A_2$. Assume $A_{\mathrm{side}}$ is the vector field interpolation contribution of the two edges:

$$\begin{gathered} A_{\mathrm{side}}={\mathit{N}}_1{A}_1+{\mathit{N}}_2{A}_2 \\ =({\mathit{N}}_1+{\mathit{N}}_2)A_1 \end{gathered}$$

(47)

Make a line integration of $A$ along this side and name it as $A_{\mathrm{unif}}$. Note that the following integration is orthogonal:

$${{\displaystyle \int }}_{eKL}{\mathit{N}}_{IL}\mathrm{d}s={\mathsf{\delta }}_{\left(IL\right)\left(KL\right)}$$

(48)

Here, $\mathrm{d}s$ is the differential vector of the line element along edge $eKL,$ and ${\delta }_{\left(IL\right)\left(KL\right)}$ represents the Kronecker delta. Thus, the integration $A_{\mathrm{unif}}$ is written as:

$$A_{\mathrm{unif}}=\underset{\mathrm{edge}1}{{\displaystyle \int }}A·\mathrm{d}s+\underset{\mathrm{edge}2}{{\displaystyle \int }}A·\mathrm{d}s=A_1+A_2=2A_1$$

(49)

Then Equation (47) is rewritten as:

$$A_{\mathrm{side}}=({\mathit{N}}_1+{\mathit{N}}_2)A_1=\left({\mathit{N}}_1+{\mathit{N}}_2\right)\frac{A_{\mathrm{unif}}}{2}={\mathit{N}}_{\mathrm{unif}}{A}_{\mathrm{unif}}$$

(50)

In Equation (50), we obtain a new shape function $N_{\mathrm{unif}}=\left(N_1+N_2\right)/2$ for the edge, as shown in Figure 5b. According to this, $N_{\mathrm{unif}}$ can be written in detail as follows:

For the serendipity quadratic hexahedron element, the edges on the sides along the $\xi$-axis direction:

$$N_k=\frac{1}{2}\left(1+4{\zeta }_k\zeta \right)\left(1-{\eta }^2\right)\left(1-{\xi }^2\right)\mathsf{\Delta }\zeta$$

(51)

The edges on the sides along $\mathrm{the} \eta$-axis direction:

$${\mathit{N}}_j=\frac{1}{8}\left(1+{\xi }_j\xi \right)\left(1+{\zeta }_j\zeta \right)\left({\xi }_j\xi +{\zeta }_j\zeta -1\right)\mathsf{\Delta }\eta$$

(52)

The edges on the sides along the $\zeta$-axis direction:

$${\mathit{N}}_k=\frac{1}{8}\left(1+{\xi }_k\xi \right)\left(1+{\eta }_k\eta \right)\left({\xi }_k\xi +{\eta }_k\eta -1\right)\mathsf{\Delta }\zeta$$

(53)

For the Lagrange quadratic hexahedron element, the edges on the sides along the $\xi$-axis direction:

$${\mathit{N}}_i=\frac{1}{8}\eta \zeta \left({\eta }_i+\eta \right)\left({\zeta }_i+\zeta \right)\mathsf{\Delta }\xi$$

(54)

The edges on the sides along the $\eta$-axis direction:

$${\mathit{N}}_j=\frac{1}{8}\xi \zeta \left({\xi }_j+\xi \right)\left({\zeta }_j+\zeta \right)\mathsf{\Delta }\eta$$

(55)

The edges on the sides along the $\zeta$-axis direction:

$${\mathit{N}}_k=\frac{1}{8}\xi \eta \left({\xi }_k+\xi \right)\left({\eta }_k+\eta \right)\mathsf{\Delta }\zeta$$

(56)

The shape functions of the edges in the element and the faces remained intact. Therefore, we obtained two new quadratic edge elements (Figure 6). There is only one edge on the sides of the two new quadratic edge elements (Figure 6); therefore, the redundant variable can be effectively eliminated.

Table 1. Elements and the distribution of their basic functions.

Element	Edges	Faces	Interior Bubbles	Total
1st	12	0	0	12
SC	24	12	0	36
SN	12	12	0	24
LC	24	24	6	54
LN	12	24	6	42

lists various types of elements and the distribution of their basis functions.

4.5. Interpolation and Convergence

Interpolation capacity and convergence tests are performed using two new types of second-order element. Those functions are approximated in the domain defined by $\Omega ={\left[-1 1\right]}^3$. Figure 7 shows the error of interpolation of the tested functions and their respective convergence rates. The two basis functions (SC, LC) have similar convergence rates, both close to one.

5. Model Studies

All MCSEM response simulations were performed on a platform with hardware and operation, as follows. The CPU was a double Intel processor, each processor having 6 cores and 12 threads, and RAM memory was 64 GB at a frequency of approximately 1333 MHz. The operating system was a 64-bit Windows 7 professional version.

5.1. Verification of the New Algorithm

To verify the new algorithm, an isotropic horizontal layered geoelectrical model, as shown in Figure 8, was considered first. In this model, air and seawater were $1000 \mathrm{m}$ thick. A $100 \mathrm{m}$ thin-layer reservoir was embedded in the sediment at $1500 \mathrm{m}$ to $1600 \mathrm{m}$ depth. The EM field response was based on the model presented in

Table 2. The conductivities of an isotropic horizontally layered geoelectric model.

Layer	$\mathbf{Conductivity} (\mathbf{s}/\mathbf{m})$
Air	${10}^{-6}$
Sea Water	$3.30$
Sediment	$1.00$
Reservoir	$0.01$

. An x-oriented $1 \mathrm{Hz}$ frequency current source was placed at the location $\left(0, 0, 950 \mathrm{m}\right)$. The receiving stations were located on the seafloor $\left(z=1000 \mathrm{m}\right)$. The modeling domain was $3 \mathrm{km} \times 3 \mathrm{km}$ square and 4 km in the vertical direction $\left(-1 \mathrm{km}\le \mathrm{z}\le 3 \mathrm{km}\right)$. A non-uniform rectangular grid was used to discretize $\mathsf{\Omega }$. We refined the grid near the current source, receiving stations, and target layer. The analytical response was calculated using the numerical integration of the Hankel transform [25,26,27,28,29].

Figure 9 and Figure 10 show the difference in the $x$-component of the electric field by applying conventional serendipity (SC)-and Lagrange (LC)-type second-order hexahedral elements, and the proposed elements (new serendipity (SN) and Lagrange element (LN)). These new quadratic edge elements were at the same level of accuracy as conventional quadratic elements, both higher than that of linear edge elements.

Table 3. The comparison of CPU time by using different methods.

Element	Number of Elements	DOF	CPU Time $\left(\mathbf{s}\right)$
1st order	332,112	867,730	30.5
SC	332,112	3,952,952	51.7
SN	332,112	2,960,408	76.2
LC	332,112	5,920,816	111.4
LN	332,112	4,924,272	221.9

lists the CPU times for different elements. Less time was required when using new quadratic elements than when using conventional quadratic elements. Therefore, one can say that the computation is reduced at no cost to accuracy.

5.2. Off-Shore Hydrocarbon ReservoirMmodel

Consider a geoelectric model with a hydrocarbon reservoir embedded in sediment, as shown in Figure 11 and Figure 12. In this model, air and seawater were $1000 \mathrm{m}$ thick, similar to the previous model. The conductivity values of each component are listed in

Table 4. The conductivities of a hydrocarbon reservoir geoelectrical model.

Layer	$\mathbf{Conductivity} \mathbf{in} \mathbf{Horizontal} (\mathbf{s}/\mathbf{m})$	$\mathbf{Conductivity} \mathbf{in} \mathbf{Vertical} (\mathbf{s}/\mathbf{m})$
Air	${10}^{-6}$	${10}^{-6}$
Sea Water	$3.30$	$3.30$
Sediment	$1.00$	$0.80$
Reservoir	$0.05$	$0.05$

. The modeling domain was a $4 \mathrm{km} \times 4 \mathrm{km}$ × 4 km cube with $1 \mathrm{km}\le z\le 3 \mathrm{km} \mathrm{vertically}$. The hydrocarbon reservoir was buried $2.0–2.1 \mathrm{km}$ deep and occupied a $1 \mathrm{km} \times 1 \mathrm{km}$ area horizontally. An $x$-oriented 1 Hz frequency current source was placed at the location $\left(-3000, 0, 950\right) \mathrm{m}$. Receiving stations were placed on the seafloor $\left(z=1000 \mathrm{m}\right)$. We used the same type of rectangular grid used in the previous model to discretize the model domain. We also refined the local grid for the current source, receiving stations, and the target layer area.

Figure 13 and Figure 14 show the $x$-component of the secondary electric field calculated with different edge elements. The results of the four types of second-order elements were in agreement with each other. Furthermore, a strong anomalous field distortion caused by the background conductivity anisotropy can also be observed.

Table 5. The cost-time of CPU when using different methods.

Element	Number of Elements	DOF	CPU Time $\left(\mathbf{s}\right)$
SC	560,000	6,679,300	141.7
SN	560,000	4,999,050	156.2
LC	560,000	9,998,100	311.4
LN	560,000	8,317,850	471.9

shows the cost-time of calculation with four types of quadratic edge elements. Similar to the previous model, the efficiencies of the two types of quadratic edge elements (SN and LN) were higher than those of conventional serendipity or Lagrange-type quadratic hexahedron elements (SC, LC).

5.3. Pyramid Type Bathymetry Model

Further verification of the new algorithm was performed by considering a pyramid-type bathymetry model (Figure 15) without a reservoir. The top area of the pyramid-type bathymetry model was $0.5 \mathrm{km}\times 0.5 \mathrm{km}$ in square, at $990 \mathrm{m}$ depth. The bottom area was $1.0 \mathrm{km}\times 1.0 \mathrm{km}$ square, at $1000 \mathrm{m}$ in depth. Air and sea water were $1000 \mathrm{m}$ thick. In the domain of pyramid-type bathymetry, the depth of seawater ranges from $900 \mathrm{m}$ to $1000 \mathrm{m}$. Figure 16 shows the pyramid-type bathymetry model with a reservoir. The domain of the reservoir was $\mathsf{\Omega }=\left\{-1000 \mathrm{m}\le x,y\le 1000 \mathrm{m},1500 \mathrm{m}\le z\le 1600 \mathrm{m}\right\}$. The modeling domain was $\mathsf{\Omega }=\left\{-4 \mathrm{km}\le x,y\le 4 \mathrm{km},-1 \mathrm{km}\le z\le 3 \mathrm{km}\right\}$. In these models, the conductivities of air, seawater, and sediment were the same as previously, whereas the reservoir’s conductivity was $0.05 \mathrm{s}/\mathrm{m}$. The electric dipole source was $x$-oriented at $\left(-3000, 0, 800\right) \mathrm{m}$ at $1 \mathrm{Hz}$ frequency, and the receiving station was located at $z=800 \mathrm{m}$.

Figure 17, Figure 18, Figure 19 and Figure 20 show the secondary electric component in the $x$-direction of the bathymetry model, which has no reservoir or one reservoir in the isotropic sediment. Figure 21, Figure 22, Figure 23 and Figure 24 show a comparison of the $x$-component of the secondary electric field of the bathymetry model with and without a reservoir in anisotropic sediment. Bathymetry and anisotropic background clearly affected the distribution of the marine controlled-source electromagnetic field (Figure 18, Figure 19, Figure 20, Figure 21, Figure 22, Figure 23 and Figure 24), which should be considered in the processing of MCSEM data. Figure 25 and Figure 26 show a comparison of the degrees of freedom (DOF) and CPU time using different elements, respectively. Clearly, these two new elements contributed significantly to reducing the CPU time in computing.

6. Conclusions

We developed a new vector finite element algorithm to solve the 3D marine CSEM modeling problem in an anisotropic medium. An algorithm based on a secondary field formulation and two new quadratic edge elements was proposed. Compared with the vector FEM algorithm based on the conventional quadratic edge element, the proposed algorithm significantly reduced the computational cost without losing accuracy. The algorithm in this paper has good generality and is suitable for solving other electromagnetic induction problems in isotropic and anisotropic media, including airborne and borehole electromagnetic methods.