Skip to content
BY 4.0 license Open Access Published by De Gruyter Open Access April 5, 2022

3D marine controlled-source electromagnetic modeling using an edge-based finite element method with a block Krylov iterative solver

Hanbo Chen, Bin Xiong, Yang Yang, Yu Han and Ziyu Cheng
From the journal Open Geosciences

Abstract

An edge-based finite element method for the numerical modeling of electromagnetic fields in complex media is presented. We used the analytical solution on an electric field in a homogeneous half space to develop a source correct factor to reduce the influence of source singularity and boundary conditions on the numerical accuracy, so that we can minimize the time required to construct the field source term in the scattered field formula. The modeling domain was discretized using an unstructured tetrahedral mesh. We transformed the complex equations of the electrical field into two real-valued equations by decomposing the field into real and imaginary components. Thereafter, we adopted a block conjugate orthogonal conjugate gradient (BL_COCR) iterative solver with an incomplete LU decomposition preconditioner, which was robust for ill-conditioned systems and efficient for multiple source electromagnetic modeling to solve the real-valued equation systems. Using the analytical solution on an electric field in a homogeneous layer model, we evaluated the accuracy of our numerical forward solution and the results showed that the source correct factor can reduce forward modeling errors associated with boundary effects and source singularities. We also applied the developed algorithm to compute the CSEM responses for typical 3D marine geo-electric models with different number of sources and compared with different iterative solvers, and the results showed that the BL_COCR solver has high computational efficiency when solving multiple right-hand term problems.

1 Introduction

The marine controlled-source electromagnetic (MCSEM) method is a well-developed tool for offshore hydrocarbon exploration [1,2]. As the depth of exploration targets and acquisition of high-quality exploration data has increased, the demand for high-accuracy MCSEM data interpretation methods have also increased [3]. Its inversion is based on forward modeling technology; therefore, research regarding highly accurate and efficient MCSEM forward modeling technology is crucial for effective inversion [4].

There are some popular numerical approaches for MCSEM forward modeling: the finite difference method (FDM) [5,6], integral equation method (IEM) [7,8], finite volume method (FVM) [9,10], and finite element method (FEM) [4,1012]. The FDM is the most straightforward and easiest to implement [5]. However, FDM formulation is generally based on rectangular grids; therefore, it cannot simulate the complex bathymetry and geo-electricity geometry well. The IEM only needs to discretize the abnormal region, but it generates a fully populated matrix. This indicates that the storage requirement and computation time will increase with an increasing number of model grid points. The FEM has the ability to simulate complex geo-electricity geometry and bathymetry because it can utilize unstructured grids to discretize the modeling domain. The FVM can effectively reduce the numerical error caused by mesh generation by the quadrature in the control volume. It can also use unstructured grids to discretize the modeling domain, but the resulting mesh generation is more difficult than that for the FEM. In this study, we used the edge-based FEM based on an unstructured mesh for 3D MCSEM modeling in an anisotropic medium.

For 3D MCSEM modeling, there are two principal approaches to solve the quasi-static variant of Maxwell’s equation. In the first method, the secondary field formulation and primary field must be computed in advance. However, to construct the field source term of the secondary field, it is necessary to calculate the primary electric field or magnetic field at nearly every discrete grid sampling point. For 3D problems, the number of electric or magnetic field sampling points in the discrete grid is more than that of the grid by a factor of approximately three and may be in the order of millions or even tens of millions. Consequently, the temporal requirements to construct a scattering field source are substantial. Increasing the number of transmitters may exceed the time required to solve the linear equation. The second method is total field formulation. Although it is not required to compute the primary field, this method requires local grid refinement of the source location to reduce the effect of source singularity, which will result in increased computation time. Therefore, to reduce the effect of sources without increasing the computation time, we applied a source correction technology to reduce the influence of the singularity of the field source and boundary conditions on the numerical accuracy [13]. Source correction technology has been applied to 2.5D direct current resistivity modeling, which has achieved positive results [13].

Subsequent to FEM analysis, the linear equation system is obtained, which may be solved using direct or iterative solvers. Direct solvers may easily solve linear equation systems with multiple right-hand sides, which is common for a typical MCSEM survey. However, it requires greater memory consumption than iterative solvers. Recently, open-source codes for controlled-source electromagnetic (CSEM) 3D forward modeling based on the FVM or FEM have been presented [1417]. Nearly all of these codes use direct solvers, such as MUMPSs or PARDISO. Using the direct method to solve the system of equations can obtain high-precision solutions, but it needs a large memory space. Compared to the direct method, the iterative solver can be easily implemented in common operating systems. However, because of the ill condition and null space property of the double curl matrix equation for the electrical field, it can be difficult to solve large linear matrix equations using traditional iterative solvers. The complex equations of the electrical field can be transformed into two real-valued equations by the decomposition of the electrical field into real and imaginary components. Subsequently, the iterative solver can solve the two resulting real-valued equations. However, for the forward modeling problem of MCSEM with multiple parameters, the traditional iterative solver operates separately on the linear system; consequently, this consumes considerable computing time. With the development of the block method, some block iterative algorithms have been proposed that solve the linear equation system with multiple right-hand sides, such as the block bi-conjugate gradient (BL_BICG), block BICG stabilized (BL_BICGSTAB) [18], block quasi-minimal residual (BL_QMR) [19], block IDR(s) (BL_IDR(s)) [20], and block GMRES (BL_GMRES) [21] methods. In this study, we introduce a new block conjugate orthogonal conjugate gradient (BL_COCR) with the residual orthonormalization technique to solve the linear equation system of multiple sources. Thus, we can expand the application of the iterative solver in marine CSEM forward modeling.

The article is divided as follows: Section 2 describes the total field formulation for 3D MCSEM modeling in anisotropic media. Boundary conditions and source term improvements are also presented in Section 2. The framework of the new block iterative Krylov subspace methods (BL_COCR) is presented in Section 2.2. We present some numerical examples that illustrate the performance of the proposed algorithm in Section 3. Finally, Section 4 presents the study conclusion.

2 Methods

2.1 Total electromagnetic field equation

The low-frequency electromagnetic field considered in geophysical applications satisfies Maxwell’s equations:

(1) × E = i ω μ 0 H ,

(2) × H = J + σ ¯ E .

where E denotes the electrical field strength, H denotes the magnetic field strength, ω denotes the angular frequency, and μ 0 denotes the magnetic permeability of free space. J represents the source current, whereas σ ¯ denotes the media conductivity tensor and is defined as follows:

(3) σ ¯ = σ x 0 0 0 σ y 0 0 0 σ z ,

where σ x , σ y , and σ z represent the spindle conductivities in the x ,   y , and  z directions, respectively.

In the practical application of MCSEM, an electric dipole source is typically used to excite the electromagnetic field. For simplicity, we only use electrical sources to simulate the response of the MCSEM. The density of source current J may be defined as follows:

(4) J = J s u δ ( x x s ) ,

where J s denotes the magnitude of the source current density, u denotes the source direction, x s denotes the source location vector, and δ denotes the Dirac delta function.

After taking the curl of equations (1) and (2), an electromagnetic wave equation formulated in terms of the electric field can be obtained as follows:

(5) × × E i ω μ 0 σ ¯ E = i ω μ 0 J .

2.2 Edge element finite element analysis

A tetrahedral mesh was used to discretize the modeling domain. Figure 1 shows a tetrahedral element with node and edge indexing. Considering a particular edge j with length l j e and with two nodes (j1 and j2) at either ends, the lowest-order edge elements are used here. The vector basis function for each edge can be defined as [22]:

(6) N j e = ( L j 1 e L j 2 e L j 2 e L j 1 e ) l j e ,

where L j 1 e and L j 2 e are the scalar basis functions for the j 1 and j 2 nodes, respectively.

Figure 1 
                  The tetrahedral element with node and edge indexing. The 1-4 is the number of the nodes. The ①-⑥ is the number of the edges.

Figure 1

The tetrahedral element with node and edge indexing. The 1-4 is the number of the nodes. The ①-⑥ is the number of the edges.

The approximated electric field corresponding to the ith edge is denoted as E i e . For a tetrahedral element, the approximated electric field E in each element can be expressed as follows:

(7) E = j = 1 6 N j E j e .

By substituting equation (7) into equation (5), and using the Galerkin’s method, the weak form of the original differential equation is obtained as follows:

(8) Ω R j = Ω N j [ × × E i ω μ 0 σ ¯ E i ω μ 0 J ] d v ,

where Ω denotes the modeling domain.

Replacing the vector basis function for the weight functions in equation (8) and expanding the equation over the entire tetrahedral mesh gives:

(9) j N edge E j e Ω ( × N i ) ( × N j ) + i ω μ 0 j N edge E j e Ω N i σ ¯ N j = i ω μ 0 j N edge Ω N i J .

After splitting the unknown coefficient, E , into real and imaginary parts, the following system is obtained:

(10) C ω μ 0 D ω μ 0 D C E R E I = 0 ω μ 0 S 1 .

where C i , j = Ω ( × N i ) ( × N j ) , D i , j = Ω N i σ ¯ N j , S 1 = Ω N j J s , E R , and E I are the real and imaginary components of the electric field E .

Before solving equation (10), a proper boundary condition must be applied. It was assumed that the boundaries of the domain are located far from the transmitter at a distance where the electromagnetic fields are negligible. Homogeneous Dirichlet boundary conditions were imposed on the outer boundary:

(11) e | Γ = 0

where Γ denotes the outer boundary domain.

2.3 Boundary conditions and source term improvement

In the CSEM, owing to the existence of an external field source, the electromagnetic field to be characterized showed a significant gradient near the sources. The singularities of the sources greatly reduce the accuracy of the numerical results. At present, there are two principal methods for reducing the influence of the field source singularity. In the first method, the scattered field formula is used and the electromagnetic field to be solved is decomposed into the superposition of the primary field, generated by the reference model (generally homogeneous space or 1D layered model), excited by both the external field source and the secondary field generated by the 3D model. Because the external source has been included in the calculation of the primary field, the secondary field source no longer has a singularity. The second method for reducing singularities is based on the total field formula of the electromagnetic field. The influence of field source singularity is reduced by grid refinement near the source location. It is worth mentioning that many authors have been working on tailored mesh design using total electric field formulation [14,15,23].

In the scattering field formula, to avoid the source singularity when calculating the secondary field, it is necessary to ensure that the resistivity of the selected reference model is consistent with that near the source. For MCSEM measurements, because the source is generally located in seawater, it is generally easy to satisfy this condition. However, for land CSEM measurements, CSEM exploration usually consists of multiple transmitter positions within a large survey area. When the excitation source is moved, if the conductivity of the region where the new excitation source is located changes, it is necessary to solve the primary electromagnetic field response of a new different background model. With the increase in the number of transmitting sources, the time consumed in the construction of the field source term increases rapidly and even exceeds the time consumed in the coefficient matrix decomposition. In contrast, in the total field formula, although the mesh refinement strategy can reduce the influence of field source singularity, it will significantly increase the number of unknowns to be solved, thus, increasing the computational overhead and reducing the solution efficiency. In addition, to satisfy the boundary conditions, whether the scattered field formula or total field formula is used, it is generally necessary to locate the outer boundary of the calculation area sufficiently far from the area of interest. The result is that the equations to be solved are very large.

To reduce the influence of field source singularity and boundary conditions on numerical accuracy, the source correction technology [13] is used to correct the source term in equation (11) of the total field control equation. This avoids the increase in calculation time caused by the local mesh refined near the sources. The correction is performed as follows:

  1. A homogenous half-space model with a conductivity of σ H that is the best initial guess of the average conductivity in the model space is first defined.

  2. The FEM algorithm described above is used to discrete the homogenous half space model to obtain the following equation:

    (12) A ( σ H ) E H = b H ,

    where ( σ H ) denotes the forward operator and b H denotes the right-hand terms of a homogeneous model and can be computed using equation (9). However, equation (11) does not need to be solved.

  3. It is also assumed that E 0 denotes the electric field of the homogenous half-space model, which can be computed analytically for the half-space background conductivity ( σ H ) [24]. If there are no errors arising from the mesh or source singularity, E 0 should satisfy the following equation:

    (13) A ( σ H ) E 0 = b 0

  4. Because of the discretization error influence and that of the boundary condition, the electric field E H will differ from the analytical solution E 0 . Therefore, source correction can be obtained using the following equation:

(14) S corr = b 0 b H = A ( σ H ) E 0 b H ,

By substituting the source correction, S corr , into equation (10), the following equation is obtained:

(15) A ( σ T ) E T = b T + S corr = A ( σ H ) E 0

where b T denotes the right-hand terms of an inhomogeneous model with a conductivity of σ T . After solving equation (15), the electric field, E T , can be obtained, which must be solved. E T is the electrical field within the model space, and σ T is the true conductivity. The modified vector b T + S corr is no longer a vector with large gradients at the source locations. It now contains additional information that corrects for the inadequacies of the forward operator A ( σ H ) at the source locations and boundaries. The corrected source term greatly reduces errors at source locations and errors associated with the boundary conditions for the solution to the true conductivity field, σ T , with minimal additional computational cost.

2.4 BL_COCR iterative solver

Multiple sources are often used to excite artificial electromagnetic fields to improve the resolution of MCSEM exploration. Therefore, in the data interpretation, the forward equations with multiple right-hand terms (multiple RHSs) must be solved. Generally speaking, there are two methods for solving such multi-source equations: direct and iterative solvers. However, when solving multi-source problems, conventional iterative solvers can only solve the forward system of each source separately, which consumes considerable computing time. With the development of the iterative Krylov solver, some block iterative Krylov subspace methods are particularly designed for solving efficient linear systems with multiple RHSs [25,26]. These require one or more matrix product operations of the form AV , with V C n × p , an arbitrary rectangular matrix, per iteration step. Thus, they can solve typical memory bottlenecks of direct methods. To the best of our knowledge, there is no literature available regarding solving the forward equation of marine CSEM with multiple sources using a block iterative solver. In this article, we introduce a new block Krylov subspace iterative solver (BL_COCR) that can be considered as an extension of the COCR algorithm proposed in previous studies [27] for simplicity.

The framework for the BL_COCR solver is shown in Table 1, which also shows that the BL_COCR solver must compute the matrix products AP m and AR m + 1 at each iteration step. Fortunately, AP m can be computed using explicit matrix multiplication. Moreover, AR m + 1 can be used to compute the recurrence relation at line 9 in Table 1; therefore, the computational complexity can be reduced significantly. To further improve the computational efficiency, the BL_COCR iterative solver was applied with an incomplete LU decomposition preconditioner to solve the linear equation system of 3D MCSEM modeling with multiple sources.

Table 1

Framework for the block Krylov subspace iterative solver

Algorithm for the BL_COCR method
1. X 0 C n × p is an initial guess, R 0 = B A X 0
2. Set P 0 = R 0 , U 0 = V 0 = A R 0
3. for m = 0,1,…, until R m F / R 0 F ε do
4. solve ( U m T U m ) α m = R m T V m for α m
5. X m + 1 = X m + P m α m
6. R m + 1 = R m U m α m and V m + 1 = A R m + 1
7. Solve ( R m T V m ) β m = R m + 1 T V m + 1 for β m
8. P m + 1 = R m + 1 + P m β m
9. U m + 1 = V m + 1 + U m β m
10. End for

3 Results

To validate our algorithm, we developed the 3D MCSEM forward modeling code using MATLAB language. The code can be run on the MATLAB 2019R platform. Although we only tested the code on MATLAB 2019R, it can also run on new version of MATLAB. The code may be run for MCSEM forward modeling without the installation of any components. It is important to note that we adopted simple adding-on methods for matrix assembly. In equation (9), we loop through parameters i and j from 1 to N d ( N d is the number of dofs in an element) for all elements and find the global indices of i and j . Non-zero entries are added to the global stiffness matrix A stored in the coordinate format. This method of integrating the global matrix is more rapid than the conventional method of cyclic calculation of the local element matrix. We used COMSOL 5.3 software to generate unstructured grids. The MCSEM responses of the examples in the study were calculated using the same platform with the following parameters: 2 Intel(R) Xeon(R) CPU E5-2630 2.40 GHz, 64GB of memory, and 64-bit Windows 7.

3.1 1D model with a semi-analytical solution

To verify the accuracy of the 3D CSEM forward algorithm developed in this study, we first considered a 1D layered model with a quasi-analytical solution. Figure 2 is a standard model used to simulate high resistivity oil and gas reservoirs in an MCSEM [2]. The conductivity of the reservoir was parameterized as 0.001 S/m and its thickness was 100 m. The top of the reservoir was located 1 km from the seafloor. The resistivity of the sedimentary layer was 1 S/m. The thickness of the seawater layer was 1 km, and its conductivity was 3.3 S/m. The source was an x-oriented electric dipole at (0, 0, 0.9) km with a frequency of 0.25 Hz. The source dipole moment was 1 Am. The receivers were placed on the seafloor. The analytical solution was computed using the Hankel transform method. The modeling domain was set as Ω = { 10  km x , y 10  km , 4  km z 8  km } . The modeling domain is discrete using a tetrahedral mesh, as shown in Figure 3. The grid was refined near the target layer and observation surface. The number of elements and edges were 263,314 and 313,523, respectively.

Figure 2 
                  Sketch of a marine-based 1D layered model.

Figure 2

Sketch of a marine-based 1D layered model.

Figure 3 
                  Tetrahedral mesh at y = 0 m for a marine-based 1D layered model. The color scale indicates the tetrahedral element size.

Figure 3

Tetrahedral mesh at y = 0 m for a marine-based 1D layered model. The color scale indicates the tetrahedral element size.

Figure 4 shows a comparison between the FEM solution and quasi-analytic solutions for the marine 1D layer model. The corrected FEM solution is in agreement with the quasi-analytical solutions. Figure 5 shows the error of the FEM solution relative to the analytical solution. When the source correction was not performed, near the transmitting dipole (the transceiver distance being less than 1 km), the electromagnetic field Ex component had a large error. The maximum error of amplitude could reach approximately 10%, and the maximum error of the phase could reach 2°, which were primarily caused by the singularity of the source. It should be noted that to reduce the size of the linear equations to be solved and reduce the amount of calculation for the iterative solution, we did not refine the mesh near the source to reduce the influence of the source singularity. When the source correction technology is adopted, the influence of source singularity upon numerical accuracy rapidly reduces. The maximum error of the amplitude is reduced to less than 5%, and the maximum error of the phase is reduced to less than 1°. In addition, with the increase in transceiver distance, the influence of source singularity will be reduced rapidly. Similarly, the error will also be reduced rapidly. The amplitude error of the electromagnetic field component will be reduced to less than 2%, and the phase error will be within 1°. Therefore, we can conclude that the source correction technology proposed in this study is feasible and effective.

Figure 4 
                  Accuracy comparison between the FEM numerical solution and quasi-analytic solution for a marine 1D layer model. (a) Amplitudes of the x component of the electric field; (b) corresponding phase.

Figure 4

Accuracy comparison between the FEM numerical solution and quasi-analytic solution for a marine 1D layer model. (a) Amplitudes of the x component of the electric field; (b) corresponding phase.

Figure 5 
                  Error comparison between the FEM numerical solution and quasi-analytic solution for a marine 1D layer model. (a) Relative errors in amplitude; (b) corresponding phase difference.

Figure 5

Error comparison between the FEM numerical solution and quasi-analytic solution for a marine 1D layer model. (a) Relative errors in amplitude; (b) corresponding phase difference.

3.2 3D marine reservoir model with a flat seafloor

In this section, we consider a model with a reservoir in anisotropic media, as shown in Figure 6. A 2 km × 2 km × 0.2 km reservoir was embedded in the sediment with a central point at (1, 0, 1.6) km. The depth of the seawater was 1,000 m. The conductivities of the air and sea water were 10 6 and 3.3 S/m, respectively. The conductivity of the reservoir was 0.05 S/m. The horizontal and vertical conductivities of the sediment were 1 and 0.8 S/m, respectively. The excitation source was an x-oriented horizontal electric dipole with a moment of 10 Am and a frequency of 1.0 Hz. It was placed at (−3, 0, 0.9) km. The receivers were placed at the seafloor (z = 1 km). The modeling domain was selected as Ω = { 20 km x , y 20 km , 4 km z 8 km } . Likewise, we used the tetrahedral mesh to discretize the modeling domain, as shown in Figure 7. To illustrate the effect of anisotropic conductivity in the data, we performed simulations for the isotropic model (with a sediment conductivity equal to 1 S/m). We compared the BL_COCR method against other block Krylov subspace methods (for example, BL_QMR, BL_BICGSTAB, and BL_IDR(4)). We used the parameters of the iterations and the final true relative residual norm, which is denoted as TRR. TRR can be obtained by the following equation: TRR = log 10 ( b A X 2 / B 2 ) . The initial guess of the iterative solution was set as X 0 = 0 C n × p , where n is the number of solutions, p is the number of RHSs, the maximum number of iterations is 8,000, and the convergence error is set to 10 11 .

Figure 6 
                  The model with a sediment reservoir. The red and green rectangles denote the electric dipole source and reservoir, respectively.

Figure 6

The model with a sediment reservoir. The red and green rectangles denote the electric dipole source and reservoir, respectively.

Figure 7 
                  The tetrahedral mesh at y = 0 for the reservoir model. The color scale indicates the tetrahedral element size.

Figure 7

The tetrahedral mesh at y = 0 for the reservoir model. The color scale indicates the tetrahedral element size.

Figure 8 shows the responses of both anisotropic and isotropic models in the form of the absolute values of the fields (top row) normalized by the corresponding results of the reference model (bottom row). The reference model is a three-layer model (air-seawater-sediment). Its geometric parameters and conductivity values are the same as those of the isotropic model. A strong distortion of the anomalous field can be observed owing to the anisotropy of the background conductivity. Figure 9 shows the convergence plot of the BL_COCR solver for 3D reservoir models (anisotropic and isotropic models) for the data with a frequency of 1 Hz compared to the convergence plots of the BL_QMR, BL_BICGSTAB, and BL_IDR(4) solvers. It may be observed that the BL_COCR solver is more stable than the BL_QMR, BL_BICGSTAB, and BL_IDR(4) solvers. Table 2 shows a comparison of the computation time of different solvers, which shows that although it requires more iterations, the BL_COCR solver is more competitive than the BL_QMR method in terms of CPU time and TRR. BL_COCR is more robust than BL_BICGSTAB and BL_IDR(4) solvers in terms of iterations, computation time, and TRR.

Figure 8 
                  Responses of both anisotropic and isotropic models in the form of the absolute values of the fields (a and b) normalized by the corresponding results of the reference model (c and d).

Figure 8

Responses of both anisotropic and isotropic models in the form of the absolute values of the fields (a and b) normalized by the corresponding results of the reference model (c and d).

Figure 9 
                  Convergence curves for different solvers for the 3D marine reservoir model with isotropic and anisotropic sediments. (a) Isotropic model; (b) anisotropic model.

Figure 9

Convergence curves for different solvers for the 3D marine reservoir model with isotropic and anisotropic sediments. (a) Isotropic model; (b) anisotropic model.

Table 2

Comparison of the numerical results of different iterative solvers for the complex salt model for a single source

Method Iteration TRR CPU time (s) Number of RHSs Models
BL_COCR 1,226 −10.13 703.8956 1 Isotropic
BL_QMR 718 −10.05 1710.8114 1 Isotropic
BL_BICGSTAB 1,388 −10.06 1362.2743 1 Isotropic
BL_IDR(4) 1,677 −10.01 1350.2598 1 Isotropic
BL_COCR 1,435 −10.15 795.7586 1 Anisotropic
BL_QMR 723 −10.07 1764.6317 1 Anisotropic
BL_BICGSTAB 1,803 −10.13 1213.9026 1 Anisotropic
BL_IDR(4) 1,940 −10.18 1399.3758 1 Anisotropic

3.3 Complex salt model

In this section, we consider a complex model based on a scaled Society of Exploration Geophysicists (SEG) salt dome [28,29]. Figure 10 shows the salt dome, which extends from −1.5 to 2 km in x and y directions and from 1.3 to 4.7 km in the vertical direction. The complex bathymetry of the complex model is shown in Figure 11. The depth of the seawater was 0.82–1.1 km. The average seawater depth is 1 km. We considered four layers as the background for the complex model: air, seawater, sediment, and basement. The conductivities of air, seawater, and sediment were the same as those of previous models. The depth of the sediment was 5 km. The horizontal and vertical conductivities of the basement were 0.01 and 0.008 S/m, respectively. The horizontal and vertical conductivities of the salt were 0.05 and 0.005 S/m, respectively.

Figure 10 
                  Scaled SEG salt dome model.

Figure 10

Scaled SEG salt dome model.

Figure 11 
                  The 3D salt home model bathymetry.

Figure 11

The 3D salt home model bathymetry.

We used seven horizontal electric dipole sources for this complex model to demonstrate the effectiveness of the BL_COCR solver for the problem with multiple right-hand sides. The sources were placed at y = 0 m and z = 0.6 km with an interval of 1 km. The frequency of the electric dipole source was 1 Hz. The receivers were placed at z = 0.8 km. The modeling domain was Ω = { 6 km x , y 6 km , 1 km z 9 km } and was discretized into a 3D unstructured mesh with 2,815,495 elements. To obtain a high accuracy FEM solution, we refined the mesh near the source and receiver locations. This refining was also applied to the bathymetry profile.

Figure 12 shows the x component of the total electric field and y component of the total magnetic field for the seven electric dipole sources. The different block algorithms were compared in terms of the number of iterations and CPU time. Figure 13 shows the convergence plot of the BL_COCR solver for the complex model for the data with a frequency of 1 Hz compared to the convergence plots of the BL_QMR, BL_BICGSTAB, and BL_IDR(4) solvers. When solving multiple right-hand term problems, BL_COCR was more stable than BL_QMR, BL_BICGSTAB, and BL_IDR(4) solvers. Table 3 shows the numerical results of different iterative solvers for the complex salt model for seven sources. Block iterative solvers can not only solve the linear equation system with multiple RHSs but may reduce the computation time when compared to that of traditional iterative solvers. In addition, the BL_COCR iterative solver computation time is less than that of some block solvers (BL_QMR, BL_BICGSTAB, and BL_IDR). Therefore, block iterative solvers can be an alternative method to solve the multiple source problems that are common in a typical CSEM survey.

Figure 12 
                  Total field for all seven sources.

Figure 12

Total field for all seven sources.

Figure 13 
                  Convergence curves for different solvers with multiple RHSs for the scaled SEG salt dome model.

Figure 13

Convergence curves for different solvers with multiple RHSs for the scaled SEG salt dome model.

Table 3

Comparison of the numerical results of different iterative solvers for the complex salt model for seven sources

Method Iteration TRR CPU time (s) Number of RHSs
ILU-COCR 7,469 −10.03 3254.8963 1
ILU-QMR 5,033 −10.01 9260.1258 1
ILU-BICGSTAB 8,543 −10.02 6169.7615 s 1
ILU-IDR(4) 9,261 −10.05 5017.7254 1
ILU-BL_COCR 3,367 −10.03 1352.9377 7
ILU-BL_QMR 2,719 −10.02 5556.0754 7
ILU-BL_BICGSTAB 4,449 −10.02 2280.7842 7
ILU-BL_IDR(4) 5,323 −10.01 2007.0925 7

4 Discussion

We developed an edge-based finite element algorithm for 3D MCSEM forward modeling based on the total electric field formulation and unstructured mesh. The tetrahedral element mesh is used to discretize the model region to ensure that geometry complex models can be simulated. Compared with the forward algorithm based on the scattered field formulation, the construction of the field source term for total electric field formulation does not need to solve the electromagnetic field of the homogeneous half space or layered model for many times. The increase in the number of the sources only causes a small increase in the amount of calculation, which is also conducive to the CSEM measurement of multiple sources. In addition, the developed correction technology can effectively improve the accuracy of the numerical solution. The error of the amplitude of electric field and phase in the region near the source can be reduced by 50%. So as to avoid the increase in calculation caused by the local mesh refined near the source. We solved the linear equation system using the BL_COCR iterative solver with incomplete LU decomposition results as the preconditioner. We found that the BL_COCR iterative solver is more efficient than conventional COCR Krylov subspace methods for the multiple-source EM problem. We also compared the efficiency of the BL_COCR iterative solver with other block solvers (e.g., BL_QMR, BL_BICGSTAB, and BL_IDR). The results showed that the BL_COCR solver is faster than those block solvers and it may be a competitive algorithm for solving linear systems with multiple right-hand sides. However, we solved the matrix equation of the real and imaginary parts of the electric field instead of the complex electric field matrix equation, which will increase the amount of calculation and memory requirements. Therefore, in the next step, we will study how to use the block iterative solver to directly solve the matrix equation satisfied by the complex electric field, so as to reduce the computational consumption. The proposed algorithm is widely applicable, has expansibility, and may be extended to other EM forward modeling problems. In the future, we will implement 3D MCSEM inversion based on the proposed forward modeling algorithm.

Abbreviations

BL_BICG

block bi-conjugate gradient

BL_BICGSTAB

block bi-conjugate gradient stabilized

BL_COCR

block conjugate orthogonal conjugate gradient

BL_GMRES

block GMRES

BL_IDR(s)

block IDR(s)

BL_QMR

block quasi-minimal residual

COCR

conjugate orthogonal conjugate gradient

CSEM

controlled-source electromagnetic modeling

EM

electromagnetic

FD

finite difference method

FEM

finite element method

FVM

finite volume method

MCSEM

marine controlled-sources electromagnetic method

RHS

right-hand terms

SEG

society of exploration geophysicists

TRR

true relative residual norm

E e

approximated electric field in each element

E i e

approximated electric field corresponding to the ith edge

ω

electric field

E T

electric field

E

electrical field strength

A ( σ H )

forward operator

σ H

half-space background conductivity

l 1 and l 2

indices of the first and second nodes for the ith edge, respectively

l i e

length of ith edge

( L 1 e , L 2 e , L 3 e , L 4 e )

linear node basis functions

H

magnetic field strength

μ 0

magnetic permeability of free space

σ ¯

media conductivity tensor

Ω

modeling domain

p

number of right-hand sides

Γ

outer boundary domain

E R and E I

real and imaginary components of the electric field, E

J

source current

J s

source current magnitude

u

source direction

x s

source location vector

σ x , σ y , σ z

spindle conductivities of the x ,   y ,   z directions, respectively

σ T

true conductivity

C = Ω ( × N i e ) ( × N j e ) d v

stiffness matrix

D = Ω e N i e σ ¯ N j e d v

mass matrix

S 1 = Ω e N i J d v

source term

b 0

right-hand terms of a homogeneous model generated by FEM

b H

right-hand terms of a homogeneous model

b T

right-hand terms of an inhomogeneous model

Ex

x component of the electric field

Acknowledgements

The authors would like to thank the reviewers and editors for valuable comments and suggestions.

  1. Funding information: This research was supported by the National Natural Science Foundation of China (No. 42174080), Natural Science Foundation of Guangxi (No. 2020gxnsfaa297079), and Research Start-up Foundation of Guilin University of Technology (No. RD2100002165), Guangxi Science and Technology Plan Innovation Team Project (GXNSFGA380004), and China Postdoctoral Science Foundation (2021MD703820).

  2. Author contributions: Conceptualization: Hanbo Chen and Bin Xiong; data curation: Ziyu Cheng; formal analysis: Hanbo Chen and Ziyu Cheng; methodology: Hanbo Chen; supervision: Bin Xiong; validation: Hanbo Chen and Yang Xiang; writing-original draft: Hanbo Chen and Yu Han; writing review and editing: Hanbo Chen, Yang Yang, Yu Han, and Ziyu Cheng. All authors have read and agreed to the published version of the manuscript.

  3. Conflict of interest: The authors state no conflict of interest.

References

[1] Srnka LJ, Carazzone JJ, Ephron MS, Eriksen EA. Remote reservoir resistivity mapping. Lead Edge. 2006;25(8):972–5.10.1190/1.2148193Search in Google Scholar

[2] Constable S, Srnka LJ. An introduction to marine controlled-source electromagnetic methods for hydrocarbon exploration. Geophysics. 2007;72(2):WA3–12.10.1190/1.2432483Search in Google Scholar

[3] Streich R. 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: direct solution and optimization for high accuracy. Geophysics. 2009;75:F95F105.10.1190/1.3196241Search in Google Scholar

[4] Cai HZ, Xiong B, Zhdanov M. Three-dimensional marine controlled-source electromagnetic modelling in anisotropic medium using finite element method. Chin J Geophys. 2015;58(8):2839–50.Search in Google Scholar

[5] Newman GA, Alumbaugh DL. Frequency‐domain modelling of airborne electromagnetic responses using staggered finite differences1. Geophys Prospect. 1995;43(8):1021–42.10.1111/j.1365-2478.1995.tb00294.xSearch in Google Scholar

[6] Davydycheva S, Druskin V, Habashy T. An efficient finite-difference scheme for electromagnetic logging in 3D anisotropic inhomogeneous media. Geophysics. 2003;68(5):1525–36.10.1190/1.1620626Search in Google Scholar

[7] Wannamaker PE, Hohmann GW, SanFilipo WA. Electromagnetic modeling of three-dimensional bodies in layered earths using integral equations. Geophysics. 1984;49(1):60–74.10.2172/5349441Search in Google Scholar

[8] Zhdanov MS, Lee SK, Yoshioka K. Integral equation method for 3D modeling of electromagnetic fields in complex structures with inhomogeneous background conductivity. Geophysics. 2006;71:333–45.10.1190/1.2358403Search in Google Scholar

[9] Haber E, Ascher UM, Aruliah DA, Oldenburg DW. Fast simulation of 3D electromagnetic problems using potentials. J Comput Phys. 2000;163(1):150–71.10.1006/jcph.2000.6545Search in Google Scholar

[10] Haber E, Ruthotto L. A multiscale finite volume method for Maxwell’s equations at low frequencies. Geophys J Int. 2014;199(2):1268–77.10.1093/gji/ggu268Search in Google Scholar

[11] Schwarzbach C, Börner RU, Spitzer K. Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics – a marine CSEM example. Geophys J Int. 2011;187(1):63–74.10.1111/j.1365-246X.2011.05127.xSearch in Google Scholar

[12] Puzyrev V, Koldan J, de la Puente J, Houzeaux G, Vázquez M, Cela JM. A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling. Geophys J Int. 2013;193(2):678–93.10.1093/gji/ggt027Search in Google Scholar

[13] Pidlisecky A, Knight R. FW2_5D: A MATLAB 2.5-D electrical resistivity modeling code. Comput Geosci. 2008;34(12):1645–54.10.1016/j.cageo.2008.04.001Search in Google Scholar

[14] Castillo-Reyes O, de la Puente J, Cela JM. PETGEM: a parallel code for 3D CSEM forward modeling using edge finite elements. Comput Geosci. 2018;119:123–36.10.1016/j.cageo.2018.07.005Search in Google Scholar

[15] Castillo-Reyes O, de la Puente J, García-Castillo LE, Cela JM. Parallel 3-D marine controlled-source electromagnetic modelling using high-order tetrahedral Nédélec elements. Geophys J Int. 2019;219(1):39–65.10.1093/gji/ggz285Search in Google Scholar

[16] Werthmüller D, Rochlitz R, Castillo-Reyes O, Heagy L. Towards an open-source landscape for 3-D CSEM modelling. Geophys J Int. 2021;227(1):644–59.10.1093/gji/ggab238Search in Google Scholar

[17] Rochlitz R, Skibbe N, Günther T. custEM: customizable finite-element simulation of complex controlled-source electromagnetic data. Geophysics. 2019;84(2):F17–33.10.1190/geo2018-0208.1Search in Google Scholar

[18] El Guennouni A, Jbilou K, Sadok H. A block version of BiCGSTAB for linear systems with multiple right-hand sides. Electron Trans Numer Anal. 2003;16(129–42):2.Search in Google Scholar

[19] Freund RW, Malhotra M. A block QMR algorithm for non-Hermitian linear systems with multiple right-hand sides. Algebra Appl. 1997;254(1–3):119–57.10.1016/S0024-3795(96)00529-0Search in Google Scholar

[20] Du L, Sogabe T, Yu B, Yamamoto Y, Zhang SL. A block IDR (s) method for nonsymmetric linear systems with multiple right-hand sides. J Comput Appl Math. 2011;235(14):4095–106.10.1016/j.cam.2011.02.035Search in Google Scholar

[21] Vital B. Etude de quelques methodes de resolution de probl´emes lin´eaires de grande taille sur multiprocesseur. PhD Thesis. Rennes: Universit´e de Rennes I; 1990.Search in Google Scholar

[22] Jin JM. The finite element method in electromagnetics. New York, United States: John Wiley & Sons; 2015.Search in Google Scholar

[23] Plessix RE, Darnet M, Mulder WA. An approach for 3D multisource, multifrequency CSEM modeling. Geophysics. 2007;72(5):SM177–84.10.1190/1.2744234Search in Google Scholar

[24] Ward SH, Hohmann GW. Electromagnetic theory for geophysical applications. In: Electromagnetic methods in applied geophysics: society of exploration geophysicists. Vol. 1, Theory. Tulsa, United States: Society of Exploration Geophysicists; 1988. p. 130–311.10.1190/1.9781560802631.ch4Search in Google Scholar

[25] Zhang J, Zhao J. A novel class of block methods based on the block AA T-Lanczos bi-orthogonalization process for matrix equations. Int J Comput Math. 2013;90(2):341–59.10.1080/00207160.2012.718072Search in Google Scholar

[26] Gutknecht MH. Block Krylov space methods for linear systems with multiple right-hand sides: an introduction. In: Siddiqi AH, Duff IS, Christensen O, editors. Modern mathematical models, methods and algorithms for real world systems. New Delhi: Anamaya Publishers; 2006. p. 420–47.Search in Google Scholar

[27] Sogabe T, Zhang SLA. COCR method for solving complex symmetric linear systems. J Comput Appl Math. 2007;199(2):297–303.10.1016/j.cam.2005.07.032Search in Google Scholar

[28] Aminzadeh F, Burkhard N, Long J, Kunz T, Duclos P. Three dimensional SEG/EAEG models – an update. Lead Edge. 1996;15(2):131–4.10.1190/1.1437283Search in Google Scholar

[29] Um ES, Harris JM, Alumbaugh DL. An iterative finite element time-domain method for simulating three-dimensional electromagnetic diffusion in earth. Geophys J Int. 2012;190(2):871–86.10.1111/j.1365-246X.2012.05540.xSearch in Google Scholar

Received: 2021-11-22
Revised: 2022-02-22
Accepted: 2022-03-17
Published Online: 2022-04-05

© 2022 Hanbo Chen et al., published by De Gruyter

This work is licensed under the Creative Commons Attribution 4.0 International License.