Fictitious Point Technique Based on Finite-Difference Method for 2.5D Direct-Current Resistivity Forward Problem

: With the widespread application of the direct-current resistivity method, searching for accurate and fast-forward algorithms has become the focus of research for geophysicists and engineers. Three-dimensional forward modeling can be the best way to identify geo-electrical anomalies but are hampered by computational limitations because of the large amount of data. A practical compromise, or even alternative, is represented by 2.5D modeling characterized using a 3D source in a 2D medium. Thus, we develop a 2.5D direct-current resistivity forward modeling algorithm. The algorithm incorporates the finite-difference approximation and fictitious point technique that can improve the efficiency and accuracy of numerical simulation. Firstly, from the boundary value problem of the electric potential generated by the point source, the discrete expressions of the governing equation are derived from the finite-difference approach. The numerical solutions of the discrete electric potential are calculated after the approximate treatment of the boundary conditions with a finite-difference method based on a fictitious point scheme. Secondly, through the simulation of a homogeneous half-space model and a one-dimensional model, and compared with the analytical results, the correctness and stability of the finite-difference forward algorithm are verified. Lastly, through the numerical simulation for a two-dimensional model, 2.5D direct-current sounding responses are summarized, which can provide a qualitative interpretation of field data.


Introduction
Direct-current resistivity is a surface geophysical method that can provide essential information for investigating subsurface geological structures by injecting electric current into the Earth and measuring the corresponding voltage.It can find applications in mineral resource exploration [1][2][3], environmental and urban geological surveys [4][5][6][7], and humanitarian geophysics [8].Regular two-dimensional (2.5D) and three-dimensional (3D) surveys are now conducted.In these multidimensional cases, the interpretation of apparent resistivity data requires an accurate modeling approach.The computer resources required for 3D direct-current resistivity simulation severely restrict the practical applicability of any automated interpretation technology.The 3D direct-current resistivity forward modeling can be expected not only from the direct development of 3D algorithms but also from a reasonable compromise between the degree of required forward model complexity and the level of computer resources consumed.A practical compromise, or even alternative, is represented by 2.5D direct-current resistivity modeling characterized by the use of a 3D point source in a 2D geological structure.Numerical modeling approaches have been developed and applied for direct-current resistivity forward modeling, such as finite-difference and finite-element methods, as the approximation of partial differential equations (Poisson equation or Helmholtz equation) [9].The finite-difference method is an effective approach for solving partial differential equations by approximating second-order derivatives with a different scheme in the governing equation [10][11][12][13].Several studies have explored the efficiency and accuracy of finite-difference methods for direct-current resistivity modeling [14][15][16][17][18].However, it is not easy to calculate the exact space-domain electrical potential.The finite-element method is another numerical approach for the direct-current resistivity modeling, which solves the variational problem in integral form for the electrical potential, and then carries out a numerical integration [19][20][21][22].The finite-element methods use complex real-world information, such as topography and bathymetry, when building density models, and improve the flexibility of 2D, 2.5D, or 3D discretization domains [23].Unfortunately, high-accuracy direct-current resistivity modeling needs fine grids, which will lead to high computational costs [24].Additionally, there are some other numerical schemes for the 2.5D direct-current resistivity modeling, such as the finite-volume method [25], the element-free Galerkin method [26], and the boundary-element method [27].These methods provide a good research basis for both numerical simulations of 2D geo-electrical models and numerical simulations in related geophysical disciplines.
In this work, the finite-difference method to incorporate a fictitious point technique was developed for 2.5D direct-current resistivity modeling, which is different to the traditional finite-difference numerical method without a fictitious point scheme.The time consumption for solving the finite-difference linear equation system without the fictitious point technique would be longer than that of the finite-difference method with the fictitious point approach.In addition, the finite-difference forward algorithm incorporated with the fictitious point technique will give more accurate approximate results than the traditional finite-difference approach, especially in dealing with the 3D direct-current resistivity forward problem.
The rest of this paper is arranged as follows.The boundary value problem for the 2.5D wavenumber domain potential is introduced in Section 2. Section 3 is devoted to presenting the high-performance implementation of the finite-difference algorithm with a fictitious point technique for 2.5D direct-current resistivity modeling.In Section 4, several numerical experiments are used to confirm the accuracy and efficiency of the forward algorithm.The main discussions and conclusions are summarized in Sections 5 and 6, respectively.Appendix A contains the MATLAB code for 2.5D direct-current resistivity modeling.Appendices B and C present the mathematical expressions of analytical solutions in homogenous half-space and two-layer resistivity models, respectively.

Statement of 2.5D Direct-Current Resistivity Forward Problem
The 3D space-domain electrical potential u(x, y, z) generated by a point source located at (x s , 0, 0) on the ground surface can be described as a Poisson equation [28]: where δ is the unit Dirac (impulse) function, Q(x, y, z) denotes the constant steady-state current density specified at a point, Ω presents the computational domain, and (x s , y s , z s ) indicate the coordinates of the point source of charge injected in the x-y-z space.At the air-earth interface Γ S , shown in Figure 1, the electrical potential satisfies the following Neumann boundary condition: also applied.It can be written as where n is the outward-pointing normal vector on the truncated boundary   .
Assuming the strike of the geological structure is aligned with the y-axis, the above problem can be simplified and efficiently simulated in the Fourier domain.The forward Fourier-cosine transform is defined as follows: where  presents the wavenumber and ( )


denotes the wavenumber domain electrical potential.Applying the forward Fourier-cosine transform to Equations ( 1)-( 3), after some calculation, the Helmholtz boundary value problem for the wavenumber domain electrical potential ( ) where S  and   denote the corresponding 2D boundary of the computational domain.
The constant steady-state current density Q in ( ) xz  space can be related to the current density I injected at ( ) where A  denotes a representative area in the x-z plane related to the injection source point ( ) On the distant boundary or truncated boundary Γ ∞ , the Neumann boundary can be also applied.It can be written as where n is the outward-pointing normal vector on the truncated boundary Γ ∞ .
Assuming the strike of the geological structure is aligned with the y-axis, the above problem can be simplified and efficiently simulated in the Fourier domain.The forward Fourier-cosine transform is defined as follows: where λ presents the wavenumber and U(x, λ, z) denotes the wavenumber domain electrical potential.Applying the forward Fourier-cosine transform to Equations ( 1)-( 3), after some calculation, the Helmholtz boundary value problem for the wavenumber domain electrical potential U(x, λ, z) can be expressed as where Γ S and Γ ∞ denote the corresponding 2D boundary of the computational domain.
The constant steady-state current density Q in (x, λ, z) space can be related to the current density I injected at (x s , z s ) by where ∆A denotes a representative area in the x-z plane related to the injection source point (x s , z s ).

Methodology
To solve the wavenumber domain electrical potential in Equation ( 5), we should divide the computational domain Ω into (N − 1) × (M − 1) rectangular elements, shown in Figure 2. The mesh is designed to be rectangular with an irregular grid spacing in the x-direction and z-direction.The nodes in the x-axis are represented by i = 1, 2, 3, • • • , N, and the nodes in the z-axis are represented by j = 1, 2, 3, • • • , M. The representative area ∆A for a source point in the interior can be written as and in the limit, at the ground surface with z → 0 , To solve the wavenumber domain electrical potential in Equation ( 5), we should di-vide the computational domain  into (N − 1) × (M − 1) rectangular elements, shown in Figure 2. The mesh is designed to be rectangular with an irregular grid spacing in the xdirection and z-direction.The nodes in the x-axis are represented by for a source point in the interior can be written as ( ) ( ) and in the limit, at the ground surface with 0 z → , ( )

Discretization of the Helmholtz Governing Equation
The 2D conductivity distribution ( )  is discretized at each node (i, j) by  and the wavenumber domain numerical solution , ij U at each node (i, j) needs to be cal- culated.At any interior node (i, j) with an irregular grid spacing in the x-axis and z-axis, the finite-difference approximation of second-order derivatives in Equation ( 5) can be written as Discretization for two-dimensional geo-electric model.

Discretization of the Helmholtz Governing Equation
The 2D conductivity distribution σ(x, z) is discretized at each node (i, j) by σ i,j and the wavenumber domain numerical solution U i,j at each node (i, j) needs to be calculated.At any interior node (i, j) with an irregular grid spacing in the x-axis and z-axis, the finitedifference approximation of second-order derivatives in Equation ( 5) can be written as and .
Substituting Equations ( 6) and ( 7) into Equation ( 5), the finite-difference equation in the discretized form at any interior node can be written as follows: Using scientific notation, Equation ( 8) can be written as where C ij W is the coupling coefficient related to the node (i − 1, j) and node (i, j), C ij E is the coupling coefficient related to the node (i, j) and node (i + 1, j), C ij S is the coupling coefficient related to the node (i, j − 1) and node (i, j), C ij S is the coupling coefficient related to the node (i, j) and node (i, j + 1), and C ij P is the self-coupling coefficient at node (i, j).All these coupling coefficients are given by

Boundary Nodes Located on the Top Edge
For all top boundary nodes (i, j) It can be implemented by assuming a fictitious row of nodes in the air at j = 0, such that the wavenumber domain electrical potential U i,2 and the conductivity σ i,2 at nodes (i, 2) are reflected at the imaginary nodes (i, 0).This fictitious point technique leads to the finite-difference form of Equation (5) given by where the coupling coefficients are given by

Boundary Nodes Located on the Bottom Edge
For all bottom boundary nodes (i, j) ∂n = 0 can be written as Using the fictitious point technique, it will lead to the finite-difference form of Equation ( 5) given by where the coupling coefficients are given by

Boundary Nodes Located on the Left Edge
Using the fictitious point technique, the difference equation for boundary nodes (i, j) with j = 2, 3, • • • , M − 1 and i = 1, has the form given as follows: where the coupling coefficients are given by

Boundary Nodes Located on the Right Edge
Using the fictitious point technique, the finite-difference equation for boundary nodes (i, j) with j = 2, 3, • • • , M − 1 and i = N, has the form given as follows: where the coupling coefficients are given by

Corner Nodes with Fictitious Point Technique
For the top-left corner node (1, 1), the finite-difference equation can be written as where For the top-right corner node (N, 1), the finite-difference equation is where (∆zj) 2 , and For the bottom-left corner node (1, M), the finite-difference equation is where (∆zj−1) 2 , and For the bottom-right node (N, M), the finite-difference equation is where (∆zj−1) 2 , and

Solution of the Discrete Forward Equation System
With the governing boundary value problem of Equation ( 5) being approximated using the finite-difference method in all N × M grid nodes, the linear algebraic equations (Equations ( 9), (15), ( 20), ( 25), (30) and ( 36)-(38)) must be properly arranged into a linear system as follows: where K is a matrix originating from the coupling coefficients, U is the unknown wavenumber domain electrical potentials vector, and p is the source vector containing I/2∆A at the current electrode positions.The matrix K generated by the finite-difference algorithm based on the fictitious point technique in Equation (39) can be a sparse, symmetric, and positive-definite matrix.Figure 3 shows the non-zero element distribution for an 8 × 8 grid (just for illustration purposes).Usually, Equation (39) can be considered an ill-conditioned problem.Therefore, the linear equation system generated by the finite-difference approximation for the wavenumber domain electrical potential can be solved by Krylov-type iterative methods.Meanwhile, the appropriate pre-conditioners can also significantly improve the speed of convergence.In forward modeling, the ILU-BICGSTAB iterative method is used, which combines a bi-conjugate gradient stabilization algorithm [29,30] with an incomplete LU decomposition [31].For the numerical simulations of a simple 2D model, the discrete elements are designed as 100 × 100 in the whole computational domain and it takes 0.12 s by the ILU-BICGSTAB iterative method.Figure 4 shows the fast convergence curve and total number of the ILU-BICGSTAB iterative method.For multiple point sources, Equation (39) needs to be solved multiple times and the apparent resistivity needs to be calculated based on the corresponding observation device.For the numerical simulations of a simple 2D model, the discrete elements are designed as 100 × 100 in the whole computational domain and it takes 0.12 s by the ILU-BICGSTAB iterative method.Figure 4 shows the fast convergence curve and total number of the ILU-BICGSTAB iterative method.For multiple point sources, Equation (39) needs to be solved multiple times and the apparent resistivity needs to be calculated based on the corresponding observation device.

Calculation of Spatial-Domain Electrical Potential and Apparent Resistivity
After solving the discrete forward equation system, we can obtain the wavenumber domain electrical potentials U(x, λ, z) on each node, and then spatial-domain electrical potentials u(x, y, z) can be calculated by the Fourier inverse transform.On the profile (y = 0), the Fourier inverse transform formula can be expressed as Using numerical integration, Equation (40) can be rewritten as where N is the number of wavenumbers; r = √ x 2 + z 2 is the distance from the measured point on the profile to the source point; λ i are the discrete wavenumbers; and g i is the corre-sponding weighting coefficient.We use the optimization scheme developed by Xu et al. [32] to select the wavenumber.The wavenumber used in our 2.5D direct-current resistivity forward modeling is shown in Table 1.
For the numerical simulations of a simple 2D model, the discrete elem signed as 100 × 100 in the whole computational domain and it takes 0.12 BICGSTAB iterative method.Figure 4 shows the fast convergence curve and of the ILU-BICGSTAB iterative method.For multiple point sources, Equati to be solved multiple times and the apparent resistivity needs to be calcula the corresponding observation device., , u x y z can be calculated by the Fourier inverse transform.On = 0), the Fourier inverse transform formula can be expressed as  In direct-current resistivity exploration, a current source +I and a current sink −I are used to energize the conductive Earth, and the apparent resistivity can be defined as follows: where G is a geometric factor, which depends upon the locations of electrodes, and ∆u is the electrical potential difference.According to the above 2.5D finite-difference algorithm analysis, the complete workflow is shown in Figure 5, which can compute the approximate numerical electrical potential and apparent resistivity.

Numerical Experiments
The Lenovo Workstation P520 with an 8-Core Intel Xeon W-2145 processor including 32 GB of RAM was applied to execute our numerical simulation.The forward algorithm was developed in MATLAB (R2023a), and the subroutine is given in Appendix A.

Benchmark with Homogeneous Half-Space
The benchmark of the accuracy of the finite-difference forward algorithm began with a homogeneous half-space model.The computational size was designed as 200 m × 100 m.The ground surface was considered flat and uniform, and homogeneous with a conductivity of 0.01 S/m.The positive and negative current electrode were placed at (−20 m, 0 m) and (20 m, 0 m), respectively, and the current was set as I = 1A.For the numerical simulation, we solved the problem on a 100 × 50 grid.
The numerical spatial-domain electrical potentials for the homogeneous half-space model are shown in Figure 6.It is evident that the spatial-domain electrical potentials computed by the finite-difference method with the fictitious point technique agree with the analytical solutions given in Appendix B. The finite-difference numerical results indicate that our fictitious point technique can provide high accuracy for 2.5D direct-current resistivity modeling.

Numerical Experiments
The Lenovo Workstation P520 with an 8-Core Intel Xeon W-2145 processor including 32 GB of RAM was applied to execute our numerical simulation.The forward algorithm was developed in MATLAB (R2023a), and the subroutine is given in Appendix A.

Benchmark with Homogeneous Half-Space
The benchmark of the accuracy of the finite-difference forward algorithm began with a homogeneous half-space model.The computational size was designed as 200 m × 100 m.The ground surface was considered flat and uniform, and homogeneous with a conductivity of 0.01 S/m.The positive and negative current electrode were placed at (−20 m, 0 m) and (20 m, 0 m), respectively, and the current was set as I = 1A.For the numerical simulation, we solved the problem on a 100 × 50 grid.
The numerical spatial-domain electrical potentials for the homogeneous half-space model are shown in Figure 6.It is evident that the spatial-domain electrical potentials computed by the finite-difference method with the fictitious point technique agree with the analytical solutions given in Appendix B. The finite-difference numerical results indicate that our fictitious point technique can provide high accuracy for 2.5D direct-current resistivity modeling.The relative root-mean-square error was used to measure the overall accuracy of finite-difference scheme [33]: Error 100 per cent where N and M are the numbers of observation nodes in the x-and y-direction, resp tively, and , ij u represent the numerical and analytical solution, respectiv Based on uniform meshes in the x-direction for the homogeneous half-space model, convergence plot in the log-log scale for the relative root-mean-square error is show Figure 7.The error decreases as the grid step size becomes smaller, and numerical stab increases with respect to grid size.This further illustrates the accuracy of our finite-diff ence forward algorithm.The relative root-mean-square error was used to measure the overall accuracy of our finite-difference scheme [33]: where N and M are the numbers of observation nodes in the xand y-direction, respectively, and u i,j and u i,j represent the numerical and analytical solution, respectively.Based on uniform meshes in the x-direction for the homogeneous half-space model, the convergence plot in the log-log scale for the relative root-mean-square error is shown in Figure 7.The error decreases as the grid step size becomes smaller, and numerical stability increases with respect to grid size.This further illustrates the accuracy of our finite-difference forward algorithm.

Two-Layer Conductivity Model
A two-layer conductivity model was applied for another numerical experiment to verify the accuracy of our forward modeling code.This two-layer model included a lowconductivity medium in the first layer and a high-conductivity medium in the second layer, as shown in Figure 8.The thickness of the low-conductivity medium was designed to be 4 m with a conductivity of 0.005 S/m, while the thickness of the high-conductivity medium was designed to be 100 m with a bearing conductivity of 0.05 S/m.The computational domain was designed as 200 m × 200 m.The direct-current resistivity data were computed using the Wenner array by varying AB/2 from 1.5 m to 75 m, and the apparent resistivity curve calculated by the finite-difference scheme is displayed in Figure 9a.The finite-difference numerical results were compared to the analytical solutions given in Appendix C from Telford et al. [34].The comparison between finite-difference numerical and analytical apparent resistivities shows a relative error below 2%, shown in Figure 9, which indicates that the 2.5D forward algorithm can satisfy practical requirements in terms of simulation accuracy.The maximum error occurs near the point source and AB/2 is approximately 70 m, but their relative errors are all within 2%.Further analysis shows that the errors in the 2.5D finite-difference forward modeling of direct-current sounding mainly come from three aspects: firstly, the errors are caused by singularity problems near the point source; secondly, the errors can be caused by truncating boundary conditions and finite-difference discretization; lastly, the errors are caused by numerical integration of the discrete Fourier inverse transform.
where N and M are the numbers of observation nodes in the x-and y-dire tively, and , i j u and , i j u represent the numerical and analytical solution Based on uniform meshes in the x-direction for the homogeneous half-spa convergence plot in the log-log scale for the relative root-mean-square erro Figure 7.The error decreases as the grid step size becomes smaller, and num increases with respect to grid size.This further illustrates the accuracy of ou ence forward algorithm.

Two-Layer Conductivity Model
A two-layer conductivity model was applied for another numerical e verify the accuracy of our forward modeling code.This two-layer model in conductivity medium in the first layer and a high-conductivity medium layer, as shown in Figure 8.The thickness of the low-conductivity medium to be 4 m with a conductivity of 0.005 S/m, while the thickness of the high medium was designed to be 100 m with a bearing conductivity of 0.05 S/m.The computational domain was designed as 200 m × 200 m.The direct tivity data were computed using the Wenner array by varying AB/2 from 1 and the apparent resistivity curve calculated by the finite-difference schem in Figure 9a.The finite-difference numerical results were compared to the a tions given in Appendix C from Telford et al. [34].The comparison between ence numerical and analytical apparent resistivities shows a relative err

A Rectangular Conductive Body in Half-Space
To illustrate the accuracy and the efficiency of our finite-difference forward algorithm, the finite-element method with the structured rectangular mesh was adopted.In order to show a comparison of accuracy and efficiency, the numerical experiment was applied for a simple 2D model, which is shown in Figure 10.It was a symmetrical, rectangular, conductive body with a resistivity of 10 Ω • m inside a uniform half-space.The background resistivity of this 2D model was designed as 100 Ω • m.The width and the thickness of the rectangular body were 4 m and 2 m, respectively.The rectangular body was located at a depth of 2 m from the ground surface.

A Rectangular Conductive Body in Half-Space
To illustrate the accuracy and the efficiency of our finite-difference forward algorithm, the finite-element method with the structured rectangular mesh was adopted.In order to show a comparison of accuracy and efficiency, the numerical experiment was applied for a simple 2D model, which is shown in Figure 10.It was a symmetrical, rectangular, conductive body with a resistivity of 10 m  inside a uniform half-space.The background resistivity of this 2D model was designed as 100 m  .The width and the thickness of the rectangular body were 4 m and 2 m, respectively.The rectangular body was located at a depth of 2 m from the ground surface.

A Rectangular Conductive Body in Half-Space
To illustrate the accuracy and the efficiency of our finite-difference forward algorithm, the finite-element method with the structured rectangular mesh was adopted.In order to show a comparison of accuracy and efficiency, the numerical experiment was applied for a simple 2D model, which is shown in Figure 10.It was a symmetrical, rectangular, conductive body with a resistivity of 10 m  inside a uniform half-space.The background resistivity of this 2D model was designed as 100 m  .The width and the thickness of the rectangular body were 4 m and 2 m, respectively.The rectangular body was located at a depth of 2 m from the ground surface.The computational size was designed as 100 m × 100 m.The Wenner-alpha array was selected for the simulation and the apparent resistivity pseudo-section calculated by The computational size was designed as 100 m × 100 m.The Wenner-alpha array was selected for the simulation and the apparent resistivity pseudo-section calculated by finite-difference scheme is displayed in Figure 11a.In order to produce the closed contour map, the resistivity distribution characteristics of the low-resistivity anomalous body (or high-conductivity anomalous body) can be qualitatively distinguished, and the spatial distribution characteristics of the anomalous body can be accurately distinguished.
finite-difference scheme is displayed in Figure 11a.In order to produce the closed contour map, the resistivity distribution characteristics of the low-resistivity anomalous body (or high-conductivity anomalous body) can be qualitatively distinguished, and the spatial distribution characteristics of the anomalous body can be accurately distinguished.
Figure 11b shows the apparent resistivity pseudo-section of the Wenner-alpha array calculated by the finite-element method from Res2dmod software (Windows version 3.03.06)[35].The finite-difference results were compared to the finite-element solutions, and the absolute error is displayed in Figure 11c.The absolute errors are all within 4 m  .By comparing the finite-difference results with the fictitious point technique and the finiteelement results, we found that the accuracy of the two ways is almost the same, and the results agree well.If the Schlumberger array was selected to simulate the vertical electrical sounding data, the apparent resistivity pseudo-section calculated by finite-difference scheme is displayed in Figure 12a.As can be seen in Figure 12, the numerical apparent resistivities calculated by the finite-difference method with the fictitious point technique and the finiteelement method with Res2dmod software are in close agreement with each other.Compared to the finite-element results (Figure 12b), the maximum absolute error for the apparent resistivity is 3.4 m  (Figure 12c).Figure 11b shows the apparent resistivity pseudo-section of the Wenner-alpha array calculated by the finite-element method from Res2dmod software (Windows version 3.03.06)[35].The finite-difference results were compared to the finite-element solutions, and the absolute error is displayed in Figure 11c.The absolute errors are all within 4 Ω • m.By comparing the finite-difference results with the fictitious point technique and the finiteelement results, we found that the accuracy of the two ways is almost the same, and the results agree well.
If the Schlumberger array was selected to simulate the vertical electrical sounding data, the apparent resistivity pseudo-section calculated by finite-difference scheme is displayed in Figure 12a.As can be seen in Figure 12, the numerical apparent resistivities calculated by the finite-difference method with the fictitious point technique and the finite-element method with Res2dmod software are in close agreement with each other.Compared to the finite-element results (Figure 12b), the maximum absolute error for the apparent resistivity is 3.4 Ω • m (Figure 12c).
x FOR PEER REVIEW 16 of 24

Discussion
In the finite-difference algorithm for 2.5D direct-current resistivity forward modeling, if the fictitious point technique was not applied, the coefficient matrix K in Equation (39) must be unsymmetrical and non-positive.Figure 13 shows the nonzero element distribution of a finite-difference coefficient matrix formed without a fictitious point technique based on an 8 × 8 grid (just for illustration purposes).The time consumption for solving the linear equation system of the finite-difference method without the fictitious point technique could be longer than that of the finite-difference method with the fictitious point approach.
Comparing the efficiency of two computing strategies, the homogeneous half-space model in Section 4.1 is adopted.The computational size is designed as 200 m × 100 m.The positive and negative current electrodes are placed at (−20 m, 0 m) and (20 m, 0 m), respectively, and the current is set as I = 1 A. The time consumption of the finite-difference method with the fictitious point technique is 0.12 s.However, the time consumption of the finite-difference method without the fictitious point technique is 0.65 s.The computational cost is mainly consumed for solving the linear equation system.In addition, the maximum relative error for the electrical potential using the fictitious point technique is 0.39%, and that without the fictitious scheme is 31.56%(Figure 14).Therefore, the accuracy and efficiency of our 2.5D direct-current resistivity forward algorithm are superior to the finite-difference method without the fictitious point technique.

Discussion
In the finite-difference algorithm for 2.5D direct-current resistivity forward modeling, if the fictitious point technique was not applied, the coefficient matrix K in Equation (39) must be unsymmetrical and non-positive.Figure 13 shows the nonzero element distribution of a finite-difference coefficient matrix formed without a fictitious point technique based on an 8 × 8 grid (just for illustration purposes).The time consumption for solving the linear equation system of the finite-difference method without the fictitious point technique could be longer than that of the finite-difference method with the fictitious point approach.
Comparing the efficiency of two computing strategies, the homogeneous half-space model in Section 4.1 is adopted.The computational size is designed as 200 m × 100 m.The positive and negative current electrodes are placed at (−20 m, 0 m) and (20 m, 0 m), respectively, and the current is set as I = 1 A. The time consumption of the finite-difference method with the fictitious point technique is 0.12 s.However, the time consumption of the finite-difference method without the fictitious point technique is 0.65 s.The computational cost is mainly consumed for solving the linear equation system.In addition, the maximum relative error for the electrical potential using the fictitious point technique is 0.39%, and that without the fictitious point scheme is 31.56%(Figure 14).Therefore, the accuracy and efficiency of our 2.5D direct-current resistivity forward algorithm are superior to the finite-difference method without the fictitious point technique.To discuss the stability of the proposed finite-difference scheme, let us also compute the condition number of the matrix K in Equation (39), since it turns out that this is a critical quantity in determining how rapidly certain iterative methods converge.The 2norm condition number is defined by   To discuss the stability of the proposed finite-difference scheme, let us also compute the condition number of the matrix K in Equation (39), since it turns out that this is a critical quantity in determining how rapidly certain iterative methods converge.The 2norm condition number is defined by To discuss the stability of the proposed finite-difference scheme, let us also compute the condition number of the matrix K in Equation (39), since it turns out that this is a critical quantity in determining how rapidly certain iterative methods converge.The 2-norm condition number is defined by where λ max and λ min are the maximum and minimum of the eigenvalues, respectively.The 2-norm condition number values as a function of the number of nodes for different meshes in the 2D wavenumber domain are shown in Figure 15.The fact that the matrix becomes more ill-conditioned as we refine the grid is responsible for the slow-down in the BICGSTAB iterative method.

Conclusions
A high-accuracy 2.5D direct-current resistivity finite-difference forward algorithm based on a fictitious point technique is developed.We present the calculation formulas o this approach and provide a successful implementation.All mathematical formulas of the forward modeling algorithm are presented and implemented in MATLAB code.The nu merical results show that the proposed algorithm has high accuracy for 2.5D direct-cur rent resistivity forward modeling, which can provide a qualitative interpretation of field data.
Unlike the conventional finite-difference method without a fictitious point approach in the spatial domain, which can lead to a linear equation with an unsymmetrical and non positive matrix, our scheme leads to a linear equation with a sparse, symmetrical positive matrix.Therefore, the computational time and the memory usage spent on solving the linear system equations should be lower than that of the conventional finite-difference method.

Figure 4 .
Figure 4. Convergence curve of the ILU-BICGSTAB iterative method with tolerance = 10 −12 .3.4.Calculation of Spatial-Domain Electrical Potential and Apparent ResistivityAfter solving the discrete forward equation system, we can obtain the wavenumber domain electrical potentials ( ),, U x z 

Figure 3 .
Figure 3. Nonzero element distribution of finite-difference coefficient matrix with an 8 × 8 grid.

Figure 4 .
Figure 4. Convergence curve of the ILU-BICGSTAB iterative method with tolerance

Figure 5 .
Figure 5.The main steps of 2.5D direct-current resistivity finite-difference forward algorithm based on fictitious point technique.

Figure 5 .
Figure 5.The main steps of 2.5D direct-current resistivity finite-difference forward algorithm based on fictitious point technique.

Figure 6 .
Figure 6.Comparison of analytical and finite-difference numerical solution of the point source tentials in the homogeneous half-space model.

Figure 6 .
Figure 6.Comparison of analytical and finite-difference numerical solution of the point source potentials in the homogeneous half-space model.

Figure 7 .
Figure 7.The convergence plot in log-log scale using relative root-mean-square error in the interval [0.1, 1].

Figure 7 .
Figure 7.The convergence plot in log-log scale using relative root-mean-square error with step sizes in the interval [0.1, 1].

Mathematics 2024 , 24 Figure 9 .
Figure 9.Comparison of analytical and finite-difference solutions for the two-layer resistivity model: (a) analytical and numerical values of apparent resistivity; (b) relative error.

Figure 9 . 24 Figure 9 .
Figure 9.Comparison of analytical and finite-difference solutions for the two-layer resistivity model: (a) analytical and numerical values of apparent resistivity; (b) relative error.

Figure 10 .
Figure 10.A two-dimensional geo-electrical model with rectangular body buried in the homogeneous half-space.

Figure 10 .
Figure 10.A two-dimensional geo-electrical model with rectangular body buried in the homogeneous half-space.

Figure 11 .
Figure 11.Numerical results for rectangular body model by different finite-element methods with Wenner-alpha electrode configuration: (a) finite-difference solution based on fictitious point technique; (b) finite-element solution based on Res2dmod software; (c) absolute error.

sFigure 11 .
Figure 11.Numerical results for rectangular body model by different finite-element methods with Wenner-alpha electrode configuration: (a) finite-difference solution based on fictitious point technique; (b) finite-element solution based on Res2dmod software; (c) absolute error.

Figure 12 .
Figure 12.Numerical results for rectangular body model by different finite-element methods with Schlumberger electrode configuration: (a) finite-difference solution based on fictitious point technique; (b) finite-element solution based on Res2dmod software; (c) absolute error.

sFigure 12 .
Figure 12.Numerical results for rectangular body model by different finite-element methods with Schlumberger electrode configuration: (a) finite-difference solution based on fictitious point technique; (b) finite-element solution based on Res2dmod software; (c) absolute error.

Figure 13 .
Figure 13.Nonzero element distribution of finite-difference coefficient matrix formed without fictitious point technique based on an 8 × 8 grid.

Figure 14 .
Figure 14.Comparison of numerical solution without fictitious point technique and numerical solution with fictitious point technique of the point source potentials in the homogeneous half-space model.

Figure 13 .
Figure 13.Nonzero element distribution of finite-difference coefficient matrix formed without fictitious point technique based on an 8 × 8 grid.

Figure 13 .
Figure 13.Nonzero element distribution of finite-difference coefficient matrix formed without fictitious point technique based on an 8 × 8 grid.

Figure 14 .
Figure 14.Comparison of numerical solution without fictitious point technique and numerical solution with fictitious point technique of the point source potentials in the homogeneous half-space model.

Figure 14 .
Figure 14.Comparison of numerical solution without fictitious point technique and numerical solution with fictitious point technique of the point source potentials in the homogeneous halfspace model.


are the maximum and minimum of the eigenvalues, respectively The 2-norm condition number values as a function of the number of nodes for differen meshes in the 2D wavenumber domain are shown in Figure15.The fact that the matrix becomes more ill-conditioned as we refine the grid is responsible for the slow-down in the BICGSTAB iterative method.

Figure 15 .
Figure 15.Condition number values as a function of the number of nodes for different meshes.

Table 1 .
Discrete wavenumbers λ i and corresponding coefficients g i of inverse Fourier transform.