Three-Dimensional Wide-Band Electromagnetic Forward Modelling Using Potential Technique

: The efﬁcacy of Krylov subspace solvers is strongly dependent on the preconditioner applied to solve the large sparse linear systems of equation for electromagnetic problems. In this study, we present a three-dimensional (3-D) plane wave electromagnetic forward simulation over a broadband frequency range. The Maxwell’s equation is solved in a secondary formulation of the Lorentz gauge coupled-potential technique. A ﬁnite-volume scheme is employed for discretizing the system of equations on a structured rectilinear mesh. We employed a block incomplete lower-upper factorization (ILU) preconditioner that is suitable for our potential formulation to enhance the computing time and convergence of the systems of equation by comparing with other preconditioners. Furthermore, we observe their effect on the iterative solvers such as the quasi-minimum residual and bi-conjugate gradient stabilizer. Several applications were used to validate and test the effectiveness of our method. Our scheme shows good agreement with the analytical solution. Notably, from the marine hydrocarbon and the crustal model, the utilisation of the bi-conjugate gradient stabilizer with block ILU preconditioner is the most appropriate. Thus, our approach can be incorporated to optimize the inversion process.


Introduction
The effectiveness and accuracy of the electromagnetic (EM) responses in heterogeneous media are vital for exploration and delineation of the geological formation [1]. The significance of the electromagnetic technique in geophysical perspective is the diverse application to a range of geological problems from shallow to deep investigation resulting from the sensitivity of the occurrence of fluids, mineralogy, etc. [2].
Numerical simulation has been a typical scheme for modeling and analyzing three-dimensional electromagnetic data [3]. Additionally, an efficient inversion of the electromagnetic data depends on the forward model algorithm. There are various numerical techniques used to solve the three-dimensional electromagnetic forward problems such as the differential equation (finite difference (FD), finite element (FE), finite volume (FV)) and the Integral equation (IE) method. The integral equation (IE) solutions avoid the discretization of the whole domain because only the anomalous domain is discretized [4][5][6]. As a result, the degree of freedom is minimal. Though, the system of linear equations is compact and dense. Several researchers have improved this method [7][8][9][10]. However, for the differential equation such as the finite difference [11][12][13][14], it has been extensively applied because of the simplicity of the approach [15]. It requires large discretization of the computational domain because it consider both the background and anomalous region. It allows for complex earth modelling. The widespread use of the finite element approach is due to the comprehensive mathematical basis [16,17]. The finite element 2 of 11 scheme can be applied to both structured [18][19][20] and unstructured grids [21,22]. The governing differential equation is formulated using the weighted residual. The Finite volume can be utilised on both structured [23] and unstructured form [24]. These unstructured grids are essential for accurate modelling of complex and irregular shaped geometry such as the bathymetry of the earth surface. In addition, the capacity to refine the grid locally in the zone of interest allows an efficient discretization of complex earth model [25,26].
The discretized numerical solution of the Helmholtz equation resulted in continuity loss at low frequency for an iterative solver [27]. we can alleviate this problem using the divergence correction [27,28]. On the other hand, we can implement the potential technique. Several investigation has adopted the use of potentials (magnetic vector A and electric scalar ϕ) in 3-D forward simulation to solve specific geo-electromagnetic problems [29][30][31][32]. Reference [33] applied the A-ϕ disintegration of the electric field using the finite volume method to solve magnetotelluric problems. Reference [22] utilized the nodal based finite element approach of the A-ϕ decomposition to simulate controlled source electromagnetic (CSEM) problems. The edge-based finite element technique was used to observe the (galvanic and inductive) effect of a controlled source electromagnetic forward model by [21].
The frequently used preconditioners such as incomplete lower-upper factorization (ILU) and Cholesky preconditioners that are nontrivial to adopt are necessary to improve the solution of linear systems [34]. However, such preconditioners are not efficient to some type of linear systems of equation. Thus, an efficient preconditioner is vital for broadband (plane-wave) electromagnetic applications as a result of developing an effective forward solver needed for the inversion routine.
In this study, we present a wide-band plane wave electromagnetic modelling problem solved with finite volume method using potentials. We included the theoretical development in secondary formulation and applied the finite volume discretization. In addition, we describe the concept of the block ILU preconditioner. To verify our approach, we simulate and compared our results with the analytical solution and also further tested the algorithm on two synthetic models. Additionally, different iterative solvers and preconditioners were also conducted to observe the most efficient with the example models.

Formulation of Maxwell's Equations in the Secondary Field
Considering Maxwell's equations for geophysical application in the quasi-static frequency domain, we have where the time dependence is e −iωt , E is the electric field, µ o indicates the magnetic permeability in the free space, ω denotes the angular frequency, H represents the magnetic field, J = σE which signifies the current density, and σ implies the conductivity [35]. By applying the curl to Equation (1) and substituting Equation (2), one can obtain the second derivative equation for the total E field. In addition, we can partition the total field into the primary and secondary fields [36]. The secondary fields satisfy Here, E p and E s are the primary and secondary electric fields respectively. the (complex valued) conductivity can be expressed as σ = ∆σ + σ b with ∆σ = σ − σ b , σ b is the known background conductivity and Ω indicates the computed domain. For the Equation (3), the primary field are computed for a layered background conductivity model. We can decompose the magnetic induction (B) field and the electric field (E) as where, A is the magnetic vector potential and ϕ denotes the electric scalar potential. By implementing the formulation in terms of secondary potentials, the curl-curl Equation (3) can be expressed in terms of the primary potentials (A p , ϕ p ) and secondary potentials (A s , ϕ s ) as Using the vector identity ∇ × ∇ × A = −∇ 2 A + ∇(∇.A), and enforcing the Lorenz gauge condition ∇.A = iµ o ωσ b ϕ (it allows the uniqueness of the magnetic vector (A) potential), we can obtain We need an additional equation to solve the Equation (7). Taking the divergence of Equation (6), which is equivalent to enforcing the continuity of current density at boundaries. Equations (6) and (8) with homogenous (Dirichlet) boundary condition (A s , ϕ s ) ≡ (0, 0) are used to obtain the unknowns.

Finite Volume Analysis
In computing numerically, the partial differential equation (PDE), it is vital to discretize the model domain. Several numerical techniques can be adopted to solve the equation, but in our case, we applied the finite volume method (Figure 1) for the discretization of the models with the magnetic vector potential defined on the edge centre and electric scalar potential at the node.
where, is the magnetic vector potential and denotes the electric scalar potential. By implementing the formulation in terms of secondary potentials, the curl-curl Equation (3) can be expressed in terms of the primary potentials ( , ) and secondary potentials ( , ) Using the vector identity ∇ × ∇ × = −∇ + ∇(∇. ) , and enforcing the Lorenz gauge condition ∇. = (it allows the uniqueness of the magnetic vector ( ) potential), we can obtain We need an additional equation to solve the Equation (7). Taking the divergence of Equation (6), we get which is equivalent to enforcing the continuity of current density at boundaries. Equations (6) and (8) with homogenous (Dirichlet) boundary condition ( , ) ≡ (0, 0) are used to obtain the unknowns.

Finite Volume Analysis
In computing numerically, the partial differential equation (PDE), it is vital to discretize the model domain. Several numerical techniques can be adopted to solve the equation, but in our case, we applied the finite volume method ( Figure 1) for the discretization of the models with the magnetic vector potential defined on the edge centre and electric scalar potential at the node. The discretized governing equation of the assembled (matrix) system of equations can be written as where K denotes the operators and the conductivity terms, e represents the unknown parameters to be determined, and b represents the source term. The effectiveness of the Krylov iterative solvers regarding its convergence is dependent on the condition number of the system matrix. Since the assembled system of Equation (9) is ill conditioned due to the conductivity contrasts between the air- The discretized governing equation of the assembled (matrix) system of equations can be written as where K denotes the operators and the conductivity terms, e represents the unknown parameters to be determined, and b represents the source term. The effectiveness of the Krylov iterative solvers regarding its convergence is dependent on the condition number of the system matrix. Since the assembled system of Equation (9) is ill conditioned due to the conductivity contrasts between the air-ground interface and the use of non-uniform discretized mesh. It is essential that we can enhance the performance of the linear system by applying a preconditioner to yield.
Here, the M −1 represent the preconditioner incorporated into the linear system of equation. In this study we use the Matlab environment to implement the block ILU preconditioner and compared with other preconditioners such as the symmetric successive over-relaxation (SSOR) and symmetric Gauss Siedel (SGS).
The block ILU preconditioner is suitable for the discretization of the partial differential equation in Equation (7) due to reordering of the coefficient matrix K in block form.
The assembled matrix can be expressed in Equation (11). Where the D denotes the block diagonal matrix of the diagonal blocks, L represent the lower triangular matrix of the sub-diagonal blocks and U denotes upper triangular matrix of the diagonal blocks. The merit of the block ILU preconditioner is that it enhances the calculation of the incomplete factorization and the triangular system solution. For the algorithm details of the block ILU, we refer the reader to [34].
Numerous sets of iterative solvers can be implemented to solve a large linear system of equations. The most often used type are the [37] (quasi-minimum residual (QMR), and [38] (bi-conjugate gradient stabilizer (BiCGStab). Both iterative solver stated above are Lanczos-based methods and they utilize least storage in computing the matrix-vector product.
The calculation of the relative residual error (err) is shown in Equation (12). where E re f is an electric field from the analytical solution, and E is the solution from the our finite volume method.
And the computation of the relative residual is indicated in Equation (13). The residual tolerance is 10 −8 .

Model 1
In verifying the accuracy of our finite volume algorithm, we synthesize an electromagnetic (EM) investigation for a homogeneous 1-D model. The background domain has a resistivity value of 1 Ω m, and a 15 km-thick layer of resistivity 100 Ω m is embedded in the background half-space (as indicated in Figure 2). A period range of 1-1000 s and E x polarization are considered for the EM modelling. The extent of the mesh is 40 km in both horizontal axes while 60 km in the vertical direction. Their number of cells used is 15, 15 and 84 in the x, y, z direction, respectively. The solid dash green lines indicate the analytical (A) solution while the red circle denotes the finite volume (FV) approach.
For this model, we compared our finite volume algorithm with the analytical solution [39]. From the plot in Figure 3a,b, it shows that the apparent resistivity and phase from our finite volume scheme are in good agreement with the analytic solution. It is also observed from Figure 3c,d that the relative error trend from the numerical solution rises at the transition zone between 2 to 200 s as a result of penetrating the conductive basement region. m, and a 15 km-thick layer of resistivity 100 Ω m is embedded in the background half-space (as indicated in Figure 2). A period range of 1-1000 s and polarization are considered for the EM modelling. The extent of the mesh is 40 km in both horizontal axes while 60 km in the vertical direction. Their number of cells used is 15, 15 and 84 in the x, y, z direction, respectively. The solid dash green lines indicate the analytical (A) solution while the red circle denotes the finite volume (FV) approach. For this model, we compared our finite volume algorithm with the analytical solution [39]. From the plot in Figure 3a,b, it shows that the apparent resistivity and phase from our finite volume scheme are in good agreement with the analytic solution. It is also observed from Figure 3c,d that the relative error trend from the numerical solution rises at the transition zone between 2 to 200 s as a result of penetrating the conductive basement region.

Model 2
We consider a model consisting of two anomalies embedded in the two-layer background earth in a marine environment as indicated in Figure 4. The resistivity of the reservoir (hydrocarbon) is 100 Ω m with a size of 5000 × 5000 × 100 m 3 while the resistivity of the salt dome is 33 Ω m and its volume 3000 × 3000 × 5000 m 3 . The incident electromagnetic wave is vertically propagated at frequency of 0.3 Hz.

Model 2
We consider a model consisting of two anomalies embedded in the two-layer background earth in a marine environment as indicated in Figure 4. The resistivity of the reservoir (hydrocarbon) is 100 Ω m with a size of 5000 × 5000 × 100 m 3 while the resistivity of the salt dome is 33 Ω m and its volume 3000 × 3000 × 5000 m 3 . The incident electromagnetic wave is vertically propagated at frequency of 0.3 Hz. A comparison of our finite volume techniques using potential formulation was carried out with the finite volume method (Direct E) from [40].   A comparison of our finite volume techniques using potential formulation was carried out with the finite volume method (Direct E) from [40].  A comparison of our finite volume techniques using potential formulation was carried out with the finite volume method (Direct E) from [40].   Furthermore, BiCGStab and QMR with block ILU and other preconditioners are examined to explore the best solver for this example. Figures 7 and 8 signify the respective comparison of the convergence for BiCGStab and QMR with different preconditioners such as block ILU, SSOR, and SGS for the marine model. As revealed in the Figures 7 and 8, the block ILU preconditioner shows the best convergence rate for both solvers among the tested preconditioners. With the same preconditioner, the BiCGStab solver outperforms the QMR, e.g., the block ILU for the BiCGStab converges at 45 iterations with computing time of 6 s while 140 iterations with a run time of 18 s for the corresponding QMR. Thus, the utilization of the block ILU preconditioner with BiCGStab iterative solver is favorable to simulate the EM responses for the marine model. Furthermore, BiCGStab and QMR with block ILU and other preconditioners are examined to explore the best solver for this example. Figures 7 and 8 signify the respective comparison of the convergence for BiCGStab and QMR with different preconditioners such as block ILU, SSOR, and SGS for the marine model. As revealed in the Figures 7 and 8, the block ILU preconditioner shows the best convergence rate for both solvers among the tested preconditioners. With the same preconditioner, the BiCGStab solver outperforms the QMR, e.g., the block ILU for the BiCGStab converges at 45 iterations with computing time of 6 s while 140 iterations with a run time of 18 s for the corresponding QMR. Thus, the utilization of the block ILU preconditioner with BiCGStab iterative solver is favorable to simulate the EM responses for the marine model.   Furthermore, BiCGStab and QMR with block ILU and other preconditioners are examined to explore the best solver for this example. Figures 7 and 8 signify the respective comparison of the convergence for BiCGStab and QMR with different preconditioners such as block ILU, SSOR, and SGS for the marine model. As revealed in the Figures 7 and 8, the block ILU preconditioner shows the best convergence rate for both solvers among the tested preconditioners. With the same preconditioner, the BiCGStab solver outperforms the QMR, e.g., the block ILU for the BiCGStab converges at 45 iterations with computing time of 6 s while 140 iterations with a run time of 18 s for the corresponding QMR. Thus, the utilization of the block ILU preconditioner with BiCGStab iterative solver is favorable to simulate the EM responses for the marine model.

Model 3
In this section, we compute the EM response for a three-dimensional (3-D) crustal model. For this model, the (two) anomalous rectangular bodies are situated at the third-layer of the (five) background horizons with resistivity values of 0.8 and 80 Ω m. The size of both anomalies are 60 × 40× 4 km 3 and 20 × 40× 4 km 3 respectively (Figure 9).

Model 3
In this section, we compute the EM response for a three-dimensional (3-D) crustal model. For this model, the (two) anomalous rectangular bodies are situated at the third-layer of the (five) background horizons with resistivity values of 0.8 and 80 Ω m. The size of both anomalies are 60 × 40 × 4 km 3 and 20 × 40 × 4 km 3 respectively (Figure 9).

Model 3
In this section, we compute the EM response for a three-dimensional (3-D) crustal model. For this model, the (two) anomalous rectangular bodies are situated at the third-layer of the (five) background horizons with resistivity values of 0.8 and 80 Ω m. The size of both anomalies are 60 × 40× 4 km 3 and 20 × 40× 4 km 3 respectively (Figure 9). The computed grid domain is discretized into 60 × 53 × 28 cells in the x, y, and z-directions. Both horizontal direction is 115 km while the vertical extent is 185 km. A plane wave source at frequency of 0.001 Hz was used for the model. Figures 10 indicate the plot of the real and imaginary part of the normalized electric field for the polarization from our potential algorithm and direct E approach. The results are in good agreement with each other. Both imaginary and real parts of the reflect the horizontal changes in the model structure with sharp changes at the resistivity discontinuities. The real and imaginary part of the magnetic field for the polarization from our potential algorithm and direct formulation are shown in Figure 11, again indicating that their results are consistent. The computed grid domain is discretized into 60 × 53 × 28 cells in the x, y, and z-directions. Both horizontal direction is 115 km while the vertical extent is 185 km. A plane wave source at frequency of 0.001 Hz was used for the model. Figure 10 indicate the plot of the real and imaginary part of the normalized electric field for the E x polarization from our potential algorithm and direct E approach. The results are in good agreement with each other. Both imaginary and real parts of the E x reflect the horizontal changes in the model structure with sharp changes at the resistivity discontinuities. The real and imaginary part of the magnetic field for the E y polarization from our potential algorithm and direct E formulation are shown in Figure 11, again indicating that their results are consistent.  Additionally, we also examined the BiCGStab, and QMR solvers with various preconditioners (Block ILU, SSOR, and SGS) on this model. Moreover, as displayed in Figures 12 and 13, the Block ILU pre-conditioner works exclusively better than the other preconditioner in this example. For this model, BiCGStab with the Block ILU preconditioner outperforms the corresponding QMR, e.g., the former one converges after 220 iterations while 395 iterations for QMR. Thus, the BiCGStab iterative solver preconditioned with the Block ILU shows a better preference for this test model.  Additionally, we also examined the BiCGStab, and QMR solvers with various preconditioners (Block ILU, SSOR, and SGS) on this model. Moreover, as displayed in Figures 12 and 13, the Block ILU pre-conditioner works exclusively better than the other preconditioner in this example. For this model, BiCGStab with the Block ILU preconditioner outperforms the corresponding QMR, e.g., the former one converges after 220 iterations while 395 iterations for QMR. Thus, the BiCGStab iterative solver preconditioned with the Block ILU shows a better preference for this test model.

Conclusions
We have developed a finite volume algorithm that solves the wide-band electromagnetic (EM) problem in the 3-D isotropic media. The concept is based on formulating the Maxwell's equation into

Conclusions
We have developed a finite volume algorithm that solves the wide-band electromagnetic (EM) problem in the 3-D isotropic media. The concept is based on formulating the Maxwell's equation into the magnetic vector ( ) and electric scalar ( ) potential secondary formulation using the Lorentz

Conclusions
We have developed a finite volume algorithm that solves the wide-band electromagnetic (EM) problem in the 3-D isotropic media. The concept is based on formulating the Maxwell's equation into the magnetic vector (A) and electric scalar (ϕ) potential secondary formulation using the Lorentz gauge correction. We implement a block ILU preconditioner (and compared with other preconditioner) which helps to improve the convergence and rapidly solved the system matrix with iterative solvers such as the quasi-minimum residual (QMR), and bi-conjugate gradient stabilizer (BiCGStab).
For the model tested, the BiCGStab technique with block ILU preconditioner signifies better performance in terms of convergence rate than the other iterative scheme. The obtained results also revealed its efficacy and accuracy in comparison with other numerical methods. Future investigation will be conducted by applying the approach into the inversion routine.