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

Authors: Tao Song; Yun Liu; Yun Wang
Journal: Mathematical Problems in Engineering
 (2017)

---

## 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.


---

## Body

## 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] (1)ρ=Dρ0DT=ρ11ρ12ρ13ρ21ρ22ρ23ρ31ρ32ρ33in which ρ0=diag⁡ρx,ρy,ρz and ρx, ρy, and ρz represent the resistivity of the three principle axes,(2)D=cos⁡γcos⁡β-sin⁡γcos⁡βsin⁡βsin⁡γcos⁡α+cos⁡γsin⁡βsin⁡αcos⁡γcos⁡α-sin⁡γsin⁡βsin⁡α-cos⁡βsin⁡αsin⁡γsin⁡α-cos⁡γsin⁡βcos⁡αcos⁡γsin⁡α+sin⁡γsin⁡βcos⁡αcos⁡βcos⁡α,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(3)v=Iρ1/22π1B,where B=(r-r0)T·ρ·(r-r0), in which r0=x0,y0,z0 denotes the location of the current point source and r=x,y,z denotes the location measurement.The 3D DC resistivity modeling may be defined by the following mixed boundary-valued problem [8, 13]:(4a)∇·σ∇v=-2Iδr-r0r,r0∈Ω,(4b)∂v∂n=0r∈Γs,(4c)σ∇vn+vr-r0cos⁡r-r0∧nB=0r,r0∈Γ∞.In this,σ is the conductive tensor and σ=(ρ)-1. Here, Ω denotes the model volume, Γs 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,v, can be written as(5)v=u0+u,where u0 denotes the primary potential due to a current point source in an anisotropic half-space and u is the secondary potential due to an anomaly. The conductivity tensor of the half-space is σp (also called background conductivity tensor) and can be written as σ=σp+σ′, where σ denotes the conductivity tensor of the earth and σ′ denotes the anomalous conductivity tensor. The boundary value problem of u in the anisotropic case can be written as(6a)∇·σ∇u+∇·σ′∇u0=0r∈Ω,(6b)∂u∂n=0r∈Γs,(6c)σ∇un+σ′∇u0n=qp-qu0-qur∈Γ∞,where qp=r-r0cos⁡((r-r0)∧n)/Bp, q=r-r0cos⁡((r-r0)∧n)/B, Bp=(r-r0)T·ρp·(r-r0), and B=(r-r0)T·ρ·(r-r0).Finally, the solution for (6a)–(6c) is equivalent to minimizing the following integral:(7)Fu=∫Ω12∇u·σ∇udΩ+∫Ω∇u·σ′∇u0dΩ+12∫Γ∞qu2dΓ+∫Γ∞q-qpu0udΓ.


## 3. Finite Element Approximation

### 3.1. Mesh Design

In Figures1(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, d, so the grid in it is the same. In coarse area, the distance of the node is calculated by d×Coefficient. 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.Figure 1
The global region subdivision with coarse area. (a) 2D global mesh design; (b) 3D global mesh design.
(a)
(b)In finite element implementation, the 3D region shown in Figure2 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.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.
(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 as(8)12∫e∇u·σ∇udΩ=12ueTK1eue,in which the element matrix can be written as(9a)K1e=136VeETσE,(9b)E=a1a2a3a4b1b2b3b4c1c2c3c4.The second term in (7) can be calculated in a similar way:(10)∫e∇u·σ′·∇u0dΩ=ueTK1e′u0e,where K1e′=1/36VeETσ′E.For the third term in (7) (the boundary integral), it can be written as (11a)∫ΓΔqu2dΓ=ueTK2eue,(11b)K2e=q·S123122121120002,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 S123 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 inK2e and the area of the side.The last term in (7) is the same as the third term and is written as(12)∫Γ∞q-qpu0udΓ=ueTK2e′u0e.Substituting (8), (10), (11a), (11b), and (12) for (7) results in(13)fu=12uTKu+uTK′u0,where u and u0 are the vectors assembled by all the element vectors ue and u0e. K and K′ are the assembled matrices using (K1e,K2e) and (K1e′,K2e′). Minimization of the functional gives the linear equation(14)Ku=K′u0.The matrix was assembled with 1D compressing storage which needs less memory and the SORCG was employed to solve the linear equations (Zhou, 2001).


## 3.1. Mesh Design

In Figures1(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, d, so the grid in it is the same. In coarse area, the distance of the node is calculated by d×Coefficient. 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.Figure 1
The global region subdivision with coarse area. (a) 2D global mesh design; (b) 3D global mesh design.
(a)
(b)In finite element implementation, the 3D region shown in Figure2 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.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.
(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 as(8)12∫e∇u·σ∇udΩ=12ueTK1eue,in which the element matrix can be written as(9a)K1e=136VeETσE,(9b)E=a1a2a3a4b1b2b3b4c1c2c3c4.The second term in (7) can be calculated in a similar way:(10)∫e∇u·σ′·∇u0dΩ=ueTK1e′u0e,where K1e′=1/36VeETσ′E.For the third term in (7) (the boundary integral), it can be written as (11a)∫ΓΔqu2dΓ=ueTK2eue,(11b)K2e=q·S123122121120002,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 S123 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 inK2e and the area of the side.The last term in (7) is the same as the third term and is written as(12)∫Γ∞q-qpu0udΓ=ueTK2e′u0e.Substituting (8), (10), (11a), (11b), and (12) for (7) results in(13)fu=12uTKu+uTK′u0,where u and u0 are the vectors assembled by all the element vectors ue and u0e. K and K′ are the assembled matrices using (K1e,K2e) and (K1e′,K2e′). Minimization of the functional gives the linear equation(14)Ku=K′u0.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ρx=ρy = 0.5 Ω·m and ρz = 2.0 Ω·m.21 electrodes were deployed in bothx- and y-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 inx- and y-directions and 20 meshes in z-direction (with the same coarse meshes) and grid coefficients are shown in Table 1.Table 1
Three coarse meshes.
Meshes
Coefficients
Mesh 1
1,1
,
1,2
,
4,8
,
16,32,64,128,256,512
Mesh 2
1,1
,
1,2
,
2,4
,
8,16,32,64,128,256
Mesh 3
1,1
,
1,2
,
4,8
,
16,32,64,128,512,1024First, we use Mesh 1 to test the solution with different element splitting and the relative errors to the analytical solution are shown in Figure3.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 Figure2(c), (b) by grid shown in Figure 2(d), and (c) by grid shown in Figure 2(e).
(a)
(b)
(c)It is clear from these tests that element splitting has a significant impact on the accuracy of numerical solutions. While Figure3(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.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.
(a)
(b)
(c)From the total field simulation in Figure4, 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ρx=ρy = 0.25 Ω·m. If the background resistivity is the same as the model, the secondary potential (u) is equal to zero which will not provide a good test case. The relative error is shown in Figure 5.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.
(a)
(b)
(c)Compared to Figure4, 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

Figure6 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 Table2. 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 locationx (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.6987From Table2, 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:(15)ρa=KΔVI,in whichρa is the apparent resistivity, K is the geometry factor of an electrode array, and ΔV is the voltage at or between the potential electrodes.Figure7 shows a two-layer model with azimuthal anisotropy.Figure 7
A two-layer model with azimuthal anisotropy.In Figure7, the principal resistivity of the first layer is ρx=ρz = 100 Ω·m and ρy=10 Ω·m and the underlaying half-space is ρx=ρz = 10 Ω·m and ρy = 1 Ω·m, for both layers, α=β=γ=0°. The apparent resistivity of pole-pole array along x-direction is shown in Figure 8, compared with Li and Spitzer’s [8] solution.Figure 8
Apparent resistivity of pole-pole array alongx-direction for model shown in Figure 7.From Figure8, we can see that our results agree well with Li and Spitzer’s. Also, the results show that, for x-direction at small electrode spacings, the apparent resistivity is approximately 31.6 Ω·m, which is the geometric mean of the resistivities of the first layer, and approximately 3.16 Ω·m at large electrode spacings, which is the geometric mean of the resistivities of the underlaying half-space [8].Figure9 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 Figure9, the resistivity along the layer interface is the same (ρx=ρy = 50 Ω·m) and the resistivity perpendicular to XY plane is ρz = 200 Ω·m. “α” is the dipping angle of the TTI media and we simulated three cases with α=0°, α=45°, and α=90°. The apparent resistivity of pole-pole array for x- and y-direction is shown in Figure 10.Figure 10
Apparent resistivities alongx- and y-direction for pole-pole array with changing dipping angles.In Figure10, the apparent resistivity in x- and y-direction is the same when α=0° while it is different when α=90°. At small electrode spacings, the apparent resistivity along x-direction approaches ρm=ρx·ρy = 100 Ω·m, which is the geometric resistivity of the first layer. Resistivity along y-direction approaches 50 Ω·m while α=90°. We have got a similar conclusion by Wang [11].


### 6.2. Model with Anomaly

Figure11 shows a 2D anisotropic model embedded in an isotropic half-space.Figure 11
A 2D anisotropic model.There are 41 electrodes along they-direction and 11 lines along the x-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.Figure 12
Dipole-dipole array apparent resistivity of model shown in Figure11. 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 r=4 m. (b) Dipole spacing r=20 m.
(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 Table1 was used for this test and we used secondary field approach in this 2D FEM algorithm.Figure13 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.Figure 13
Comparison of 3D and 2D FEM numerical solutions of the model shown in Figure11. (a) Apparent resistivity pseudosections along the central line by 3D FEM. (b) Apparent resistivity pseudosections by 2D FEM.
(a)
(b)The final model is an anisotropic cube embedded in a half-space (Figure14(a)). The resistivity along the principal axis is ρx=ρy = 10 Ω·m and ρz = 100 Ω·m and the dipping angle α is in Y-Z plane. The model and survey setup are shown in Figure 14(b).Figure 14
(a) A 3D anisotropic cube in a homogeneous isotropic half-space. The principal resistivity of the cube is given byρx/ρy/ρz = 10/10/100 Ω⋅m and the background resistivity is 100 (Ω⋅m). (b) 21 survey lines were deployed along both x- and y-direction and the solid lines outline the boundary of the buried cube.
(a)
(b)Figure15 shows the pseudosection of 3D anisotropic model at lines X11 and X07 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 ρ1=10×100≈31.6 Ω·m.Figure 15
Apparent resistivity pseudosection of the model shown in Figure14. (a) Apparent resistivity pseudosection of line X11 for the isotropic model. (b) Apparent resistivity pseudosection of line X07 for the isotropic model. (c) Apparent resistivity pseudosection of line X11 for the anisotropic model. (d) Apparent resistivity pseudosection of line X07 for the anisotropic model.
(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α=45° and 90°; the pseudosections of apparent resistivity are shown in Figures 16 and 17, respectively.Figure 16
Pseudosection of the apparent resistivity for dipole-dipole array when the dipping angleα=45°. (a) Pseudosection of apparent resistivity at X11. (b) Pseudosection of apparent resistivity at X07. (c) Pseudosection of apparent resistivity at Y11. (d) Pseudosection of apparent resistivity Y07.
(a)
(b)
(c)
(d)Figure 17
Pseudosection of the apparent resistivity for dipole-dipole array when the dipping angleα=90°. (a) Pseudosection of apparent resistivity at X11. (b) Pseudosection of apparent resistivity at X07. (c) Pseudosection of apparent resistivity at Y11. (d) Pseudosection of apparent resistivity at Y07.
(a)
(b)
(c)
(d)For the case of dipping angleα=45°, the pseudosection along x-direction is not symmetrical anymore. Pseudosections at X11 and X7 are very different from that of α=0° and the shape of the low resistivity feature is indicative of the direction of dipping. While pseudosections long Y-direction are symmetrical, sections at the same distance from Y11 become quite different.Similarly, apparent resistivity sections forα=90° are very different from that of α=0° 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.


## 6.1. Layered Model

Figure6 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 Table2. 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 locationx (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.6987From Table2, 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:(15)ρa=KΔVI,in whichρa is the apparent resistivity, K is the geometry factor of an electrode array, and ΔV is the voltage at or between the potential electrodes.Figure7 shows a two-layer model with azimuthal anisotropy.Figure 7
A two-layer model with azimuthal anisotropy.In Figure7, the principal resistivity of the first layer is ρx=ρz = 100 Ω·m and ρy=10 Ω·m and the underlaying half-space is ρx=ρz = 10 Ω·m and ρy = 1 Ω·m, for both layers, α=β=γ=0°. The apparent resistivity of pole-pole array along x-direction is shown in Figure 8, compared with Li and Spitzer’s [8] solution.Figure 8
Apparent resistivity of pole-pole array alongx-direction for model shown in Figure 7.From Figure8, we can see that our results agree well with Li and Spitzer’s. Also, the results show that, for x-direction at small electrode spacings, the apparent resistivity is approximately 31.6 Ω·m, which is the geometric mean of the resistivities of the first layer, and approximately 3.16 Ω·m at large electrode spacings, which is the geometric mean of the resistivities of the underlaying half-space [8].Figure9 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 Figure9, the resistivity along the layer interface is the same (ρx=ρy = 50 Ω·m) and the resistivity perpendicular to XY plane is ρz = 200 Ω·m. “α” is the dipping angle of the TTI media and we simulated three cases with α=0°, α=45°, and α=90°. The apparent resistivity of pole-pole array for x- and y-direction is shown in Figure 10.Figure 10
Apparent resistivities alongx- and y-direction for pole-pole array with changing dipping angles.In Figure10, the apparent resistivity in x- and y-direction is the same when α=0° while it is different when α=90°. At small electrode spacings, the apparent resistivity along x-direction approaches ρm=ρx·ρy = 100 Ω·m, which is the geometric resistivity of the first layer. Resistivity along y-direction approaches 50 Ω·m while α=90°. We have got a similar conclusion by Wang [11].


## 6.2. Model with Anomaly

Figure11 shows a 2D anisotropic model embedded in an isotropic half-space.Figure 11
A 2D anisotropic model.There are 41 electrodes along they-direction and 11 lines along the x-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.Figure 12
Dipole-dipole array apparent resistivity of model shown in Figure11. 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 r=4 m. (b) Dipole spacing r=20 m.
(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 Table1 was used for this test and we used secondary field approach in this 2D FEM algorithm.Figure13 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.Figure 13
Comparison of 3D and 2D FEM numerical solutions of the model shown in Figure11. (a) Apparent resistivity pseudosections along the central line by 3D FEM. (b) Apparent resistivity pseudosections by 2D FEM.
(a)
(b)The final model is an anisotropic cube embedded in a half-space (Figure14(a)). The resistivity along the principal axis is ρx=ρy = 10 Ω·m and ρz = 100 Ω·m and the dipping angle α is in Y-Z plane. The model and survey setup are shown in Figure 14(b).Figure 14
(a) A 3D anisotropic cube in a homogeneous isotropic half-space. The principal resistivity of the cube is given byρx/ρy/ρz = 10/10/100 Ω⋅m and the background resistivity is 100 (Ω⋅m). (b) 21 survey lines were deployed along both x- and y-direction and the solid lines outline the boundary of the buried cube.
(a)
(b)Figure15 shows the pseudosection of 3D anisotropic model at lines X11 and X07 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 ρ1=10×100≈31.6 Ω·m.Figure 15
Apparent resistivity pseudosection of the model shown in Figure14. (a) Apparent resistivity pseudosection of line X11 for the isotropic model. (b) Apparent resistivity pseudosection of line X07 for the isotropic model. (c) Apparent resistivity pseudosection of line X11 for the anisotropic model. (d) Apparent resistivity pseudosection of line X07 for the anisotropic model.
(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α=45° and 90°; the pseudosections of apparent resistivity are shown in Figures 16 and 17, respectively.Figure 16
Pseudosection of the apparent resistivity for dipole-dipole array when the dipping angleα=45°. (a) Pseudosection of apparent resistivity at X11. (b) Pseudosection of apparent resistivity at X07. (c) Pseudosection of apparent resistivity at Y11. (d) Pseudosection of apparent resistivity Y07.
(a)
(b)
(c)
(d)Figure 17
Pseudosection of the apparent resistivity for dipole-dipole array when the dipping angleα=90°. (a) Pseudosection of apparent resistivity at X11. (b) Pseudosection of apparent resistivity at X07. (c) Pseudosection of apparent resistivity at Y11. (d) Pseudosection of apparent resistivity at Y07.
(a)
(b)
(c)
(d)For the case of dipping angleα=45°, the pseudosection along x-direction is not symmetrical anymore. Pseudosections at X11 and X7 are very different from that of α=0° and the shape of the low resistivity feature is indicative of the direction of dipping. While pseudosections long Y-direction are symmetrical, sections at the same distance from Y11 become quite different.Similarly, apparent resistivity sections forα=90° are very different from that of α=0° 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.



---

*Source: 8027616-2017-05-15.xml*