Research Article  Open Access
Finite Element Method for Modeling 3D Resistivity Sounding on Anisotropic Geoelectric Media
Abstract
A 3D DC finite element method forward program is developed in this paper for anisotropic geoelectric media. Both total and secondary field approaches have been implemented. In this paper, we focused on the structured grid scheme. The modeling shows that the symmetry of the structure grid determines the symmetry of the response potential around the point source for an anisotropic halfspace. Through numerical modeling with three kinds of coarse meshes by the total field approach and secondary field approach separately, a higher accuracy can be achieved via the secondary field approach. And the relative error between the numerical solution and the analytical solution is less than 2%. Modeling results contrasting with previous scholars also verify the correctness of the algorithm. Then, the numerical results of 3D models with anisotropic properties were presented and compared to models with isotropic properties. These results clearly illustrated the strong effect of anisotropy and the problems in interpretation if anisotropy was not properly addressed. This work will establish the foundation for our future effort of building a 3D inversion program with arbitrary anisotropy.
1. Introduction
The direct current (DC) method is well established under the assumption of 3D isotropic geoelectric media. Pridmore et al. [1] first discussed the application of the finite element method (FEM) approach to the 3D resistivity problem in isotropic media and a finite element method (FEM) based 3D inversion scheme was implemented by Sasaki [2]. Many techniques have been reported ([3–5]; Zhao & Yedlin 1996) for reducing computation time and improving accuracy. Recently, much effort has been devoted to the 3D problem in anisotropic media. Li and Uren [6] derived an analytical solution for point source potential in anisotropic halfspace and solutions for a layered formation with arbitrary anisotropy were also reported by Yin and Weidelt [7]. Using rectangular hexahedral elements, Li and Spitzer [8] presented a study on 3D resistivity modeling with arbitrary anisotropy. Due to the limitations of the parallelepiped element, it is challenging to model complex geoelectric structures with the approach of Li and Spitzer [8]. Zhou et al. [9] used Gaussian quadrature grid (GQG) method to model 2.5D/3D DC in heterogeneous, anisotropic media, and Greenhalgh et al. [10] developed explicit expressions for sensitivity functions in resistivity imaging of heterogeneous and fully anisotropic earth.
Recently, Wang [11] proposed an unstructured FEM to model arbitrary anisotropic resistivity and achieved good accuracy. However, the unstructured grid proposed in that paper will be very difficult to implement in an inversion program. In this paper, we describe a 3D algorithm for modeling resistivity structures with arbitrary anisotropy. We used a tetrahedral element which can be easily implemented in the inversion program and investigated the effects of grids schemes and coarse meshes on the accuracy of numerical solutions. We presented the results of a few 3D anisotropic models to investigate the influence of anisotropic media. This work will establish the foundation for our effort to develop a 3D resistivity inversion program with arbitrary anisotropy.
2. Governing Equations and Boundary Conditions
In anisotropic media, the resistivity tensor can be written as [12] in which and , , and represent the resistivity of the three principle axes,and , , and represent the three Euler angles.
According to Li and Uren [6], the electrical potential from a current point source in anisotropic halfspace can be written aswhere , in which denotes the location of the current point source and denotes the location measurement.
The 3D DC resistivity modeling may be defined by the following mixed boundaryvalued problem [8, 13]:
In this, is the conductive tensor and . Here, denotes the model volume, is the surface of the model, and is the artificial boundary where the mixed boundary condition is applied.
When the singularity removal procedure is applied to solve the secondary potential problem, the total potential, , can be written aswhere denotes the primary potential due to a current point source in an anisotropic halfspace and is the secondary potential due to an anomaly. The conductivity tensor of the halfspace is (also called background conductivity tensor) and can be written as , where denotes the conductivity tensor of the earth and denotes the anomalous conductivity tensor. The boundary value problem of in the anisotropic case can be written aswhere , , , and .
Finally, the solution for (6a)–(6c) is equivalent to minimizing the following integral:
3. Finite Element Approximation
3.1. Mesh Design
In Figures 1(a) and 1(b), the red line area denotes the computation area in which we are interested, and the dashed black line denotes the coarse area. Generally, the distance of the node in computation area is fixed, for example, , so the grid in it is the same. In coarse area, the distance of the node is calculated by . The larger the distance from the computation area to the node, the greater the coefficient, which leads to the boundary being far away from the computation area.
(a)
(b)
In finite element implementation, the 3D region shown in Figure 2 is divided into a large number of bricks first, and then with the subdivision shown in Figures 2(a) and 2(b) the bricks consist of five tetrahedral elements [2]. Figure 2(e) shows the subdivision which combined scheme A and scheme B. To obtain the symmetric grid, the mesh near the surface uses the two schemes A and B alternately, and along the vertical direction the two schemes are used alternately too.
(a)
(b)
(c)
(d)
(e)
3.2. Element Analysis
Based on Sasaki [2] and Bing and Greenhalgh [5], the difference between their deductions is that the conductivity is a tensor which is symmetric.
The first term in (7) can be written asin which the element matrix can be written as
The second term in (7) can be calculated in a similar way:where .
For the third term in (7) (the boundary integral), it can be written as in which represents the side of the tetrahedron element which coincides with the artificial boundary, the elemental matrix coefficient can be calculated as (11b), and is the area of triangle 123.
For other sides of the tetrahedron element that coincides with the boundary, the boundary integral can be calculated in a similar way but changing the corresponding elements in and the area of the side.
The last term in (7) is the same as the third term and is written as
Substituting (8), (10), (11a), (11b), and (12) for (7) results inwhere and are the vectors assembled by all the element vectors and . and are the assembled matrices using and . Minimization of the functional gives the linear equation
The matrix was assembled with 1D compressing storage which needs less memory and the SORCG was employed to solve the linear equations (Zhou, 2001).
4. Total Field Approach
To test the accuracy of the program, a Homogeneous Vertical Transverse Isotropy (VTI) halfspace model which has an analytical solution is used. The three principal resistivities of the halfspace are = 0.5 Ω·m and = 2.0 Ω·m.
21 electrodes were deployed in both  and direction with a uniform spacing of 1 meter.
We tested the total field approach first and used three different coarse meshes. There are twelve coarse meshes in  and directions and 20 meshes in direction (with the same coarse meshes) and grid coefficients are shown in Table 1.

First, we use Mesh 1 to test the solution with different element splitting and the relative errors to the analytical solution are shown in Figure 3.
(a)
(b)
(c)
It is clear from these tests that element splitting has a significant impact on the accuracy of numerical solutions. While Figure 3(c) has better accuracy, the errors produced by the total field solution are unacceptable. Using the same model, we also tested the program on three meshes shown in Table 1. The test results are shown in Figure 4.
(a)
(b)
(c)
From the total field simulation in Figure 4, it is clear that element splitting and mesh design have a significant impact on the accuracy of FEM solution and a different approach has to be investigated to achieve better accuracy.
5. Secondary Field Approach
Using the same model and same parameters, we repeated the tests using secondary field approach. To avoid a solution of zero, we set background resistivity to be = 0.25 Ω·m. If the background resistivity is the same as the model, the secondary potential () is equal to zero which will not provide a good test case. The relative error is shown in Figure 5.
(a)
(b)
(c)
Compared to Figure 4, we can see that the secondary field solution is much better in accuracy in all tested meshes and the relative error is less than 2%. In the examples below, Mesh 1 and the secondary field approach are adopted.
6. Geological Model Test
6.1. Layered Model
Figure 6 shows a twolayer Vertical Transverse Isotropy model which we used to test the accuracy of our program.
The FEM results and the analytical solutions are shown in Table 2. We use a unit current source (1 A) and the electrodes’ location represents the distance between a measuring point and a current source.

From Table 2, it is clear that FEM has achieved a very good accuracy and the relative error to the analytical solution is less than 1%. In the following case, we compute the apparent resistivity (see (15)) which is commonly used as a way to present data:
in which is the apparent resistivity, is the geometry factor of an electrode array, and is the voltage at or between the potential electrodes.
Figure 7 shows a twolayer model with azimuthal anisotropy.
In Figure 7, the principal resistivity of the first layer is = 100 Ω·m and =10 Ω·m and the underlaying halfspace is = 10 Ω·m and = 1 Ω·m, for both layers, . The apparent resistivity of polepole array along direction is shown in Figure 8, compared with Li and Spitzer’s [8] solution.
From Figure 8, we can see that our results agree well with Li and Spitzer’s. Also, the results show that, for direction at small electrode spacings, the apparent resistivity is approximately Ω·m, which is the geometric mean of the resistivities of the first layer, and approximately Ω·m at large electrode spacings, which is the geometric mean of the resistivities of the underlaying halfspace [8].
Figure 9 shows a twolayer Tilted Transverse Isotropy (TTI) model.
In Figure 9, the resistivity along the layer interface is the same ( = 50 Ω·m) and the resistivity perpendicular to plane is = 200 Ω·m. “” is the dipping angle of the TTI media and we simulated three cases with , , and . The apparent resistivity of polepole array for  and direction is shown in Figure 10.
In Figure 10, the apparent resistivity in  and direction is the same when while it is different when . At small electrode spacings, the apparent resistivity along direction approaches = 100 Ω·m, which is the geometric resistivity of the first layer. Resistivity along direction approaches Ω·m while . We have got a similar conclusion by Wang [11].
6.2. Model with Anomaly
Figure 11 shows a 2D anisotropic model embedded in an isotropic halfspace.
There are 41 electrodes along the direction and 11 lines along the direction. We used a uniform 1 m for both line and electrode spacing. To illustrate the effects of anisotropy, we simulated the responses of the same structure with/without anisotropy and the results are shown in Figure 12. For small dipoledipole spacing (4 m), the impact of anisotropy is fairly small since these short spaced electrodes sample only shallow subsurfaces.
(a)
(b)
However, for large dipole spacing (20 m), the difference can reach 50% between anisotropic and isotropic media which highlights the impact of anisotropy of the structure.
To test the accuracy of our 3D program in calculating 2D anisotropic structure, we develop a 2D version of anisotropic FEM with rectangle element and linear interpolation.
Mesh 1 shown in Table 1 was used for this test and we used secondary field approach in this 2D FEM algorithm.
Figure 13 shows the pseudosection of apparent resistivity along the central line by 3D FEM and 2D FEM. Solutions from these two programs agree very well and this test shows that accurate 2D solution can be produced from our 3D program.
(a)
(b)
The final model is an anisotropic cube embedded in a halfspace (Figure 14(a)). The resistivity along the principal axis is = 10 Ω·m and = 100 Ω·m and the dipping angle is in plane. The model and survey setup are shown in Figure 14(b).
(a)
(b)
Figure 15 shows the pseudosection of 3D anisotropic model at lines and by dipoledipole array. For comparison, the pseudosections of apparent resistivity at the same lines for an isotropic model are also shown. The resistivity of the isotropic cube is given by Ω·m.
(a)
(b)
(c)
(d)
These simulations show the significant difference in apparent resistivity between anisotropic and isotropic models and an attempt to interpret these data without considering anisotropy will certainly lead to wrong conclusions.
Simulations were also run for the case of dipping angles and 90°; the pseudosections of apparent resistivity are shown in Figures 16 and 17, respectively.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
For the case of dipping angle , the pseudosection along direction is not symmetrical anymore. Pseudosections at and are very different from that of and the shape of the low resistivity feature is indicative of the direction of dipping. While pseudosections long direction are symmetrical, sections at the same distance from become quite different.
Similarly, apparent resistivity sections for are very different from that of and become symmetrical again. While it is difficult to predict the actual dipping angle and axial resistivity from these sections, the correlations between them are very instructive. Ultimately, we will have to invert these data to resolve the target geometry and anisotropy. These modeling exercises illustrate the high sensitivity of resistivity data to anisotropy and geometry of subsurface targets and highlight the importance of anisotropy in DC resistivity data interpretation.
7. Conclusions
We implemented a 3D DC resistivity modeling program which is capable of simulating the electric field of arbitrary anisotropy. We adopted finite element method in building the system and selected symmetric tetrahedron in element splitting. To achieve better numerical accuracy, we tested various splitting and meshing schemes. By comparing our results to that of analytical solutions of simple anisotropic models, we found that symmetric grids and secondary field approach would be the best choices.
With these implementations, we tested various models with anisotropy. These results verified our program and illustrated the impact of anisotropy on DC apparent resistivity data. Without proper consideration of anisotropic effects, we would likely come to wrong conclusions in interpreting the apparent resistivity data. Finally, the forward program will serve as a solid foundation for our effort to build a 3D anisotropic inversion program.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors wish to acknowledge the workshop of GEM Chengdu 2015 (http://library.seg.org/doi/abs/10.1190/GEM2015096) while parts of this work (total field approach, secondary field approach, and model with anomaly) have been presented in this workshop. The authors also wish to acknowledge the assistance of Dr. Shenghui Li. This work was supported by the National Natural Science Foundation of China (41425017) and the National Key Research and Development Program of ChinaGeophysical Comprehensive Exploration and Information Extraction of Deep Mineral Resources (2016YFC0600505).
References
 D. F. Pridmore, G. W. Hohmann, S. H. Ward, and W. R. Sill, “An investigation of finiteelement modeling for electrical and electromagnetic data in three dimensions,” Geophysics, vol. 46, no. 7, pp. 1009–1024, 1981. View at: Publisher Site  Google Scholar
 Y. Sasaki, “3D resistivity inversion using the finiteelement method,” Geophysics, vol. 59, no. 11, pp. 1839–1848, 1994. View at: Publisher Site  Google Scholar
 T. Lowry, M. B. Allen, and P. N. SHive, “Singularity removal: a refinement of resistivity modeling techniques,” Geophysics, vol. 54, no. 3, pp. 700–707, 1989. View at: Google Scholar
 K. Spitzer, “A 3D finitedifference algorithm for DC resistivity modelling using conjugate gradient methods,” Geophysical Journal International, vol. 123, no. 3, pp. 903–914, 1995. View at: Publisher Site  Google Scholar
 Z. Bing and S. A. Greenhalgh, “Finite element threedimensional direct current resistivity modelling: accuracy and efficiency considerations,” Geophysical Journal International, vol. 145, no. 3, pp. 679–688, 2001. View at: Publisher Site  Google Scholar
 P. Li and N. F. Uren, “Analytical solution for the point source potential in an anisotropic 3D halfspace. I. Twohorizontallayer case,” Mathematical and Computer Modelling, vol. 26, no. 5, pp. 9–27, 1997. View at: Publisher Site  Google Scholar  MathSciNet
 C. Yin and P. Weidelt, “Geoelectrical fields in a layered earth with arbitrary anisotropy,” Geophysics, vol. 64, no. 2, pp. 426–434, 1999. View at: Publisher Site  Google Scholar
 Y. Li and K. Spitzer, “Finite element resistivity modelling for threedimensional structures with arbitrary anisotropy,” Physics of the Earth and Planetary Interiors, vol. 150, no. 13, pp. 15–27, 2005. View at: Publisher Site  Google Scholar
 B. Zhou, M. Greenhalgh, and S. A. Greenhalgh, “2.5D/3D resistivity modelling in anisotropic media using Gaussian quadrature grids,” Geophysical Journal International, vol. 176, no. 1, pp. 63–80, 2009. View at: Publisher Site  Google Scholar
 S. A. Greenhalgh, B. Zhou, M. Greenhalgh, L. Marescot, and T. Wiese, “Explicit expressions for the Fréchet derivatives in 3D anisotropic resistivity inversion,” Geophysics, vol. 74, no. 3, pp. F31–F43, 2009. View at: Publisher Site  Google Scholar
 W. Wang, Study of 3D arbitrary anisotropic resistivity modeling and interpretation using unstructured finite element methods [Ph.D. thesis], University of Science and Technology of China, Hefei, China, 2013.
 C. Yin, “Geoelectrical inversion for a onedimensional anisotropic model and inherent nonuniqueness,” Geophysical Journal International, vol. 140, no. 1, pp. 11–23, 2000. View at: Publisher Site  Google Scholar
 A. Dey and H. F. Morrison, “Resistivity modelling for arbitrarily shaped twodimensional structures,” Geophysical Prospecting, vol. 27, no. 1, pp. 106–136, 1979. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2017 Tao Song et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.