# Three-Dimensional Inversion of the Time-Lapse Resistivity Method on the MPI Parallel Algorithm
## 2 Rationale and Theory  
### Three-Dimensional Forward Modeling Algorithm  
In this study, the finite difference method is used to implement the three-dimensional forward modeling of the resistivity method. According to the forward modeling theory of the resistivity method derived from Zhu [11], the partial differential equation satisfied by the point source resistivity method is as follows:  
\[\
abla[\boldsymbol{\sigma}(x,y,z)\
abla\boldsymbol{\varphi}(x,y,z)]=-j\delta( x-x_{0})\delta(y-y_{0})\delta(z-z_{0}), \tag{1}\]  
where \(\boldsymbol{\sigma}\) is the electrical conductivity of rock ore (S\(\cdot\)m\({}^{-1}\)), \(\boldsymbol{\varphi}\) is the potential array (V), \(j\) is the electric current density (A\(\cdot\)m\({}^{-2}\)), and \(\delta\) is a Dirac function, that specifies the position of the supply point.  
Boundary conditions need to be added at the boundary of infinity and the interface between the surface and air in the meshed region. Boundary conditions are divided into three categories: Dirichlet boundary conditions, Neumann boundary conditions, and mixed boundary conditions.  
\[\left\{\begin{array}{l}\varphi=0\\ \sigma\frac{\partial\boldsymbol{\varphi}}{\partial n}=0\\ \frac{\partial\boldsymbol{\varphi}}{\partial n}+\frac{\cos\theta}{r}\varphi=0 \end{array}\right., \tag{2}\]  
where \(n\) is the normal vector outside the boundary, \(r\) is the distance from the source point to the boundary, and \(\theta\) is the angle between \(n\) and \(r\).  
The Dirichlet boundary condition means that the potential value on the boundary is 0. The Neumann condition boundary condition refers to the normal vector of the current density; that is, the normal derivative of the potential is 0. The mixed boundary condition not only maintains the continuity of the potential at the boundary, but also eliminates the reflection effect on the virtual source at the underground boundary [12, 13].The finite difference method is applied to discretize (1), resulting in a large sparse linear system [11], whose matrix form can be expressed as:  
\[A\boldsymbol{\varphi}=\boldsymbol{b}, \tag{3}\]  
where \(\boldsymbol{A}\) is the coefficient matrix, which is used to calculate the sparse matrix of each element through the finite difference method and then assemble it into a large coefficient matrix of all the elements; \(\boldsymbol{b}\) is the vector that contains information with respect to the field source.  
The potential is obtained by solving (3) using the BICGSTAB method. After the potential is obtained, the apparent resistivity is obtained according to (4).  
\[\left\{\begin{array}{l}\mathbf{U}=\mathbf{\alpha}\mathbf{\varphi}\\ \rho_{\mathbf{a}}=K\frac{\Delta\mathbf{U}}{l}\end{array}\right., \tag{4}\]  
where \(\mathbf{\alpha}\) is the vector that is designed to select the potential at the receiver among all grid node potentials, \(\Delta\mathbf{U}\) is the potential difference between the receiving electrodes (V), \(I\) is the electric current (A), and \(K\) is the geometric factor.  
In this study, the accuracy and precision of the finite difference method are verified. A vertical contact band model is designed, as shown in Figure 1, in which the dielectric resistivity on the left is 100 \(\Omega\cdot\)m, and the dielectric resistivity on the right is 20 \(\Omega\cdot\)m. Data acquisition is carried out using a three-pole device, where A represents the transmitting electrode, M and N represent the receiving electrode, and observations are made from left to right.  
Figure 2 shows the apparent resistivity curve of the vertical contact zone model. The numerical solution obtained by numerical simulation using the finite difference method is close to the analytical solution, and the error of both methods is less than 5% at all measurement points, which verifies the accuracy and precision of the forward numerical simulation results in this paper.