Research on numerical simulation method for 3D marine direct current resistivity

As one of the most important methods in marine electromagnetic exploration, the marine direct current (DC) resistivity method has been widely used in the exploration of marine minerals, oil and gas, as well as in large-scale engineered geological exploration and other fields. Efficient and high-precision forward algorithms are the foundation for an accurate quantitative interpretation of the marine DC resistivity method. We have implemented an efficient and high-precision forward modeling approach for 3D marine DC resistivity. First, the wave number domain anomalous potential control equations are obtained by performing a 2D Fourier transform along the horizontal direction based on the differential control equations satisfied by the marine DC anomalous potential. Instead of directly solving the 3D numerical simulation problem, we transform the 3D numerical simulation problem into multiple 1D numerical simulation problems in the wavenumber domain by dimensionality reduction for computational efficiency. Second, the boundary conditions of the governing equations are given to obtain the corresponding boundary value problems, and the anomalous potential is solved using the 1D finite element method in the wavenumber domain. Next, we perform a 2D inverse Fourier transform on the wave-number domain anomalous potential to obtain the spatial domain anomalous potential. Furthermore, compact operators are used to iteratively modify the potential and obtain high-precision numerical solutions. Finally, we demonstrate the correctness of the proposed algorithm’s solution strategy by using a hierarchical ocean model, and the efficiency of the proposed algorithm by using a spherical model.


Introduction
Marine DC resistivity methods, which use detection devices to supply direct current to seafloor formations and observe the distribution of DC electric fields in seafloor formations, achieve the goals of delineating seafloor sulfide deposits and detecting seafloor permafrost and deep oil and gas exploration.In contrast to other marine electromagnetic methods, the marine DC resistivity method can avoid the effects of electromagnetic fields in the ionosphere and near the coast.The conductivity of seawater is relatively balanced by the effects of salinity and temperature.The environment in which the electrodes are located is stable and the impedance between the electrodes and the seawater is small, which can provide a weak field source and a large current, thereby improving the measurement accuracy.In addition, the measurement of submarine cables by hauling can improve the efficiency of construction and the range of exploration.
The marine DC resistivity methods obtain the electrical distribution of subsurface structures by inversion of the underlying information observed on the seabed, and efficient and high-precision forward algorithms are a prerequisite for reliable inversion of subsurface structures.The traditional numerical simulation methods for ocean direct current resistivity mainly include integral equation method [1], finite difference method [2], finite element method [3], finite volume method [4], etc.The subsurface structure and seabed topography are extremely complex due to the influence of the sedimentary environment and mineral molecular structure.In order to achieve accurate exploration, it is necessary to finely dissect the model, resulting in a sharp increase in computational complexity.Therefore, a large number of scholars have continuously explored strategies to improve the efficiency of traditional numerical simulation algorithms.However, in general, traditional numerical simulation methods directly solve 3D problems, leading to high computational complexity and high storage requirements.To further explore and improve the computational efficiency and accuracy of numerical simulation algorithms, we have implemented a 3D numerical simulation algorithm.This algorithm uses a 2D Fourier transform to transform a 3D numerical simulation problem in the spatial domain into multiple 1D numerical problems in the wavenumber domain for solution, which enables accurate characterization and flexible dissection of the model.This approach is characterized by high efficiency and accuracy.

2.Theory
Based on the differential equation of marine direct current potential, a differential equation based on anomalous potential is derived, that is  is the anomalous conductivity, E is the total electric field vector, and  is the Laplace operator.
Perform a 2D Fourier transform along the horizontal direction on equation (1) to obtain the differential equation satisfied by the anomalous potential in the wavenumber domain, ( ) Equations ( 4) and (3) are used to form the boundary value problem of abnormal potential satisfaction in wave number domain.Discretize the solution area and use a shape function based on onedimensional quadratic interpolation to represent , , , , , a b a ax ay az U J J J  within each element.At the same time, use the 1D finite element method to solve and obtain abnormal potentials in the wavenumber domain.Perform a 2D inverse Fourier transform on the wavenumber domain anomalous potential to obtain the spatial domain anomalous potential, and add the background potential to obtain the total potential.Then, use a compact operator [5] for iterative modification to obtain high-precision numerical solutions.It should be noted that we used the standard FFT method for 2D inverse Fourier transform.
It should be noted that the selection of wave number only needs to be strictly selected according to

3.Example analysis
To verify the correctness of the above algorithm, we tested it using a two-layer horizontal model of marine with analytical solutions.The model is shown in Figure 1, where the first layer is seawater with a depth of 50 m and a resistivity of 0.3 m  .The second layer is the seabed medium with a resistivity of 1 m   .We use a pole-pole device for measurement, with the power supply electrode located 5 m above the seabed and a power supply current of 10A.The measurement electrode is arranged along the x-axis direction of the seabed.The model adopts uniform grid division, with a horizontal grid spacing of 1m and a vertical grid spacing of 0.5m.The total number of grids is 101×101×81.
The comparison results between the analytical and numerical solutions are shown in Figure 2. From it, it can be seen that the maximum relative error between the numerical and analytical solutions obtained by the proposed algorithm is 0.29%, which indicates that the high computational accuracy of the proposed algorithm.In order to verify the computational efficiency of the proposed algorithm, a spherical marine model was used for testing, and the efficiency of the finite element method [7] was compared.As shown in figure 3, there is a sphere with a radius of 10 m on the seabed, with a resistivity of 0.6 m  and a distance of 10 m from the top of the sphere to the seabed surface.The other parameters and measurements are consistent with the parameters of the previous example model.The model is divided into uniform unit grids with a horizontal step size of 2m and a vertical step size of 1m, with a total number of 201×201×101.The parameters of the personal computer used for testing were an Intel i9 9980XE with 3 GHz core and 128GB RAM.The maximum potential error obtained by the algorithm and the finite element method is 1.6%, and the maximum memory and time comparisons are shown in table 1.It can be seen that the maximum memory and time of the proposed algorithm is much smaller than that of the finite element method, indicating that the proposed algorithm is computationally efficient.

4.Conclusion
We have presented a 3D numerical simulation method for marine DC resistance, which is different from the conventional numerical simulation method.In this paper, we use the 2D Fourier transform method is used to transform the 3D numerical simulation problem into a 1D numerical simulation problem in the wavenumber domain for solution, which avoids solving large sparse matrix equations with low computational complexity and storage requirements.The algorithm has been validated by examples to be efficient and high-precision.


where b is the background conductivity, a U is the anomalous potential, a J is the scattering current density vector, which satisfies the condition of aa  = JE , b

Uk
where a represents the abnormal potential in the wavenumber domain, and , ax ay JJ and az J represent the scattering current densities in the three directions in the wavenumber domain, which satisfy = are the spatial wavenumbers in the x and y directions, respectively, with i being the imaginary unit.The boundary conditions of equation (2

Figure 1 Figure 2
Figure 1 Schematic diagram of a two-layer ocean model

Figure 3
Figure 3 Schematic diagram of a spherical marine model

Table 1
Comparison of maximum memory and time between the proposed algorithm and the finite element method