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 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 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 aswhere 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 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) 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.
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 two-layer 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 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.
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.
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.
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.
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)
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 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 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)
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/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).