Mathematical Problems in Engineering
Volume 2017 (2017), Article ID 8027616, 12 pages
https://doi.org/10.1155/2017/8027616

Research Article

Finite Element Method for Modeling 3D Resistivity Sounding on Anisotropic Geoelectric Media

Tao Song,1,2 Yun Liu,¹ and Yun Wang¹

¹State Key Laboratory of Ore Deposit Geochemistry, Institute of Geochemistry, Chinese Academy of Sciences, Guiyang 550081, China
²University of Chinese Academy of Sciences, Beijing 100049, China

Correspondence should be addressed to Tao Song; songtao@mail.gyig.ac.cn

Received 11 January 2017; Accepted 9 April 2017; Published 15 May 2017

Academic Editor: Ben T. Nohara

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.

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 half-space. 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 half-space 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 half-space can be written as

where , 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 boundary-valued 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 as

where denotes the primary potential due to a current point source in an anisotropic half-space and is the secondary potential due to an anomaly. The conductivity tensor of the half-space 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 as

where , , , 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)

Figure 1: The global region subdivision with coarse area. (a) 2D global mesh design; (b) 3D global mesh design.

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)

Figure 2: Subdivision of the region. (a) Five-tetrahedron scheme A; (b) five-tetrahedron scheme B; (c) 3D global mesh with scheme A; (d) 3D global mesh with scheme B; (e) 3D global mesh with schemes A and B.

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 as

in 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 in

where 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) half-space model which has an analytical solution is used. The three principal resistivities of the half-space 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.

Table 1: Three coarse meshes.

Meshes
Coefficients

Mesh 1

Mesh 2

Mesh 3

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)

Figure 3: Forward modeling accuracy in percentage for three kinds of subdivision relative to analytical solutions by total potential approach: (a) by grid shown in Figure 2(c), (b) by grid shown in Figure 2(d), and (c) by grid shown in Figure 2(e).

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)

Figure 4: Forward modeling accuracy in percentage for three kinds of coarse meshes relative to analytical solutions by total potential approach: (a) by Mesh 1, (b) by Mesh 2, and (c) by Mesh 3.

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)

Figure 5: Forward modeling accuracy in percentage relative to analytical solutions with three kinds of sparse meshes by secondary potential approach: (a) by Mesh 1, (b) by Mesh 2, and (c) by Mesh 3.

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 two-layer Vertical Transverse Isotropy model which we used to test the accuracy of our program.

Figure 6: A two-layer model.

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.

Table 2: Comparison of FEM and analytical results.

Electrodes location (m)
Analytical solutions (V)
3D FEM results (V)
Relative error (%)

1
0.3399
0.3430
0.9120

2
0.2538
0.2561
0.9062

3
0.2188
0.2205
0.7770

4
0.1968
0.1982
0.7012

5
0.1805
0.1818
0.6925

6
0.1678
0.1687
0.5602

8
0.1478
0.1489
0.7240

10
0.1331
0.1340
0.6987

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 two-layer model with azimuthal anisotropy.

Figure 7: A two-layer model with azimuthal anisotropy.

In Figure 7, the principal resistivity of the first layer is = 100 Ω·m and =10 Ω·m and the underlaying half-space is = 10 Ω·m and = 1 Ω·m, for both layers, . The apparent resistivity of pole-pole array along -direction is shown in Figure 8, compared with Li and Spitzer’s [8] solution.

Figure 8: Apparent resistivity of pole-pole array along -direction for model shown in Figure 7.

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 half-space [8].

Figure 9 shows a two-layer Tilted Transverse Isotropy (TTI) model.

Figure 9: A two-layer model consisting of a TTI overburden layer with variable dip over an isotropic half-space.

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 pole-pole array for - and -direction is shown in Figure 10.

Figure 10: Apparent resistivities along - and -direction for pole-pole array with changing dipping angles.

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 half-space.

Figure 11: A 2D anisotropic model.

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 dipole-dipole spacing (4 m), the impact of anisotropy is fairly small since these short spaced electrodes sample only shallow subsurfaces.

(a)

(b)

Figure 12: Dipole-dipole array apparent resistivity of model shown in Figure 11. The resistivity of the isotropic model is described as follows: the resistivity of the half-space is 20 Ω⋅m and the resistivity of the anomaly is 1 Ω⋅m. (a) Dipole spacing  m. (b) Dipole spacing  m.

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)

Figure 13: Comparison of 3D and 2D FEM numerical solutions of the model shown in Figure 11. (a) Apparent resistivity pseudosections along the central line by 3D FEM. (b) Apparent resistivity pseudosections by 2D FEM.

The final model is an anisotropic cube embedded in a half-space (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 14: (a) A 3D anisotropic cube in a homogeneous isotropic half-space. The principal resistivity of the cube is given by = 10/10/100 Ω⋅m and the background resistivity is 100 (Ω⋅m). (b) 21 survey lines were deployed along both - and -direction and the solid lines outline the boundary of the buried cube.

Figure 15 shows the pseudosection of 3D anisotropic model at lines and by dipole-dipole 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)

Figure 15: Apparent resistivity pseudosection of the model shown in Figure 14. (a) Apparent resistivity pseudosection of line for the isotropic model. (b) Apparent resistivity pseudosection of line for the isotropic model. (c) Apparent resistivity pseudosection of line for the anisotropic model. (d) Apparent resistivity pseudosection of line for the anisotropic model.

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)

Figure 16: Pseudosection of the apparent resistivity for dipole-dipole array when the dipping angle . (a) Pseudosection of apparent resistivity at . (b) Pseudosection of apparent resistivity at . (c) Pseudosection of apparent resistivity at . (d) Pseudosection of apparent resistivity .

(a)

(b)

(c)

(d)

Figure 17: Pseudosection of the apparent resistivity for dipole-dipole array when the dipping angle . (a) Pseudosection of apparent resistivity at . (b) Pseudosection of apparent resistivity at . (c) Pseudosection of apparent resistivity at . (d) Pseudosection of apparent resistivity at .

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/GEM2015-096) 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 China-Geophysical 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 finite-element modeling for electrical and electromagnetic data in three dimensions,” Geophysics, vol. 46, no. 7, pp. 1009–1024, 1981.

Y. Sasaki, “3-D resistivity inversion using the finite-element method,” Geophysics, vol. 59, no. 11, pp. 1839–1848, 1994.

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.

K. Spitzer, “A 3-D finite-difference algorithm for DC resistivity modelling using conjugate gradient methods,” Geophysical Journal International, vol. 123, no. 3, pp. 903–914, 1995.

Z. Bing and S. A. Greenhalgh, “Finite element three-dimensional direct current resistivity modelling: accuracy and efficiency considerations,” Geophysical Journal International, vol. 145, no. 3, pp. 679–688, 2001.

P. Li and N. F. Uren, “Analytical solution for the point source potential in an anisotropic 3-D half-space. I. Two-horizontal-layer case,” Mathematical and Computer Modelling, vol. 26, no. 5, pp. 9–27, 1997.

C. Yin and P. Weidelt, “Geoelectrical fields in a layered earth with arbitrary anisotropy,” Geophysics, vol. 64, no. 2, pp. 426–434, 1999.

Y. Li and K. Spitzer, “Finite element resistivity modelling for three-dimensional structures with arbitrary anisotropy,” Physics of the Earth and Planetary Interiors, vol. 150, no. 1-3, pp. 15–27, 2005.

B. Zhou, M. Greenhalgh, and S. A. Greenhalgh, “2.5-D/3-D resistivity modelling in anisotropic media using Gaussian quadrature grids,” Geophysical Journal International, vol. 176, no. 1, pp. 63–80, 2009.

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.

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 one-dimensional anisotropic model and inherent non-uniqueness,” Geophysical Journal International, vol. 140, no. 1, pp. 11–23, 2000.

A. Dey and H. F. Morrison, “Resistivity modelling for arbitrarily shaped two-dimensional structures,” Geophysical Prospecting, vol. 27, no. 1, pp. 106–136, 1979.