Abstract
The problem of electrical sounding of a medium with ground surface relief is modelled using the integral equations method. This numerical method is based on the triangulation of the computational domain, which is adapted to the shape of the relief and the measuring line. The numerical algorithm is tested by comparing the results with the known solution for horizontally layered media with two layers. Calculations are also performed to verify the fulfilment of the “reciprocity principle” for the 4-electrode installations in our numerical model. Simulations are then performed for a two-layered medium with a surface relief. The quantitative influences of the relief, the resistivity ratios of the contacting media, and the depth of the second layer on the apparent resistivity curves are established.
1. Introduction
The resistivity method is one of the oldest methods for surveying media to study the earth’s structure [1, 2]. This method is still used extensively because it provides greater depth than other methods of electrical sounding and because the mathematical methods for interpreting its results are more developed [3–12]. The method of electrical resistivity tomography (ERT), which is currently used in geophysics, is a modern modification of the resistivity method. This method is being further developed through equipment improvements, automation of the measurement process, and the possibility of conducting high density measurements [13, 14].
The common mathematical model of electrical tomography is based on Maxwell’s equations for a direct current in a conductive medium. The closing relation of this model is Ohm’s law, which states that the current density in a medium is the product of the conductivity of the medium and the gradient of the electric field potential.
In the process of electrical sounding, a pair of source electrodes introduce low-frequency currents into the ground. Then, the potential differences between the pairs of measuring electrodes are measured. Deviation from the expected pattern of potential differences for a homogeneous medium provides information about the electrical properties of the subsurface inhomogeneities. The distribution of source and measuring electrodes and the type of installation are determined by the predefined measurement protocol.
The measurement results are influenced by three main factors: the distribution of electrical conductivity in the medium, the ground surface relief, and the type of electrode array. The ultimate goal of interpreting the results is to establish the electrical conductivity distribution in the ground. This study numerically simulates the electrical sounding curves in a two-layered medium with a ground surface relief. The effects of the source electrode positions and the resistivity ratio of the contacting media on the apparent resistivity curves are calculated. The influence of the relief on the interpretation results was studied in [9–12] using finite element and finite difference methods. This paper is a continuation of our previous research on this topic [15, 16]. The novelty of our approach is that the simulations are conducted via numerical solutions of the system of integral equations. The method of integral equations is based on the theory of potential for solutions to the Laplace equation. Compared with other numerical methods, this method is highly accurate and efficient for calculating the three-dimensional potential distribution [17]. A characteristic feature of the considered models is the presence of a sharply expressed contact boundary in the medium, unlike the models used in [9–12].
The numerical results presented below are obtained for a two-layered medium with a relief in the form of a single 2D shaft. However, the results can be generalized for more complex models.
2. Mathematical Model
Our numerical models consider a medium with a piecewise constant two-dimensional distribution of electric conductivity. Let the medium comprise two layers with electric conductivities and (corresponding to specific resistivities and , resp.). Assume that the lower layer is disposed horizontally at depth and that the upper layer has an embossed surface with a single shaft (Figure 1). In media with piecewise constant conductivity, the potential of the electric field satisfies Laplace’s equation in the domains of constant conductivity and the conditions for continuity of the charge flow across boundaries between media of different resistivities.
In practice, the voltage is applied to the medium through two source electrodes, designated and . However, based on the principle of superposition of electric fields, we consider only the potential created by the point source .
Let be the ground surface and let be the contact boundary between the media with conductivities and . The potential field is the sum of the simple layer potential defined on the boundaries and the potential of point source in the top boundary :
The factor appears here as the dimensional scale for the dimensionless functions of the simple layer potential.
The integral terms in (1) correspond to the simple layer potentials of charge densities and . According to the properties of the simple layer potential, satisfies Laplace’s equation in the domain except at the borders and and at point (Figure 1). The function must satisfy the boundary condition at the inner boundary between two materials and on the surface of the medium . The first of these conditions enforces the conservation of the current through the contact boundary and is written as
Indices indicate that the normal derivatives are taken from the side of the medium with the corresponding conductivity.
Now, we can write (2) in terms of the simple layer potentials:
We consider the discontinuity in the normal derivatives of the simple layer potential on the two sides of the surface. The jump in the normal derivative of the potential taken inside and outside of the region is equal to the density of a simple layer multiplied by [19]. If the normal to the boundary is directed from medium 2 to medium 1, then
In this case, the following equalities hold:
The derivatives of the simple layer potential at the second boundary are discontinuous, but the derivatives of the potential of the other sources remain continuous. Therefore, the continuity condition for a normal current passing through a surface is written as
Substituting (5) into (6), we obtain
Hence, the following integral equation is obtained:
Here, denotes a reflection coefficient at the boundary between media 1 and 2. Substituting the potential for source electrode and into (8), we obtain
The boundary condition at the medium’s ground surface is derived from the condition that other normal currents, except for the current from the point source, do not flow into medium 1. Then, we can write the expression for the normal current that enters medium 1:
The current flowing into the medium is expressed in terms of the normal derivative of the potential, taken from the side of the medium. Due to the properties of the simple layer potential [19], the normal derivative of the potential has a discontinuity on the surface of . Let the normal be directed to the outside of medium 1. According to [19], its value on the inner side of the surface is expressed as
Substituting (11) into (10), we obtain
The above equation yields another integral equation:
Thus, we have two integral relations, (9) and (13), for the two unknown functions of the simple layer densities and , which are defined at the ground surface and at the contact boundary , respectively.
3. Numerical Method
The numerical solution of the system of integral equations is conducted by discretizing the integrals in (9) and (13). We limit the infinite surfaces and by their finite parts representing the oval-shaped regions. The boundary is mapped on the plane onto which the irregular triangular mesh is constructed. The grid is condensed along the axis at the part corresponding to the measuring line, where the field potentials are computed. Using the difference in the potentials values along the profile, the apparent resistivity of the medium is calculated, as is customary in geophysical experiments.
A typical triangulation of the calculation domain is shown in Figure 2, and the unknown functions and are computed at the triangulation nodes. The integrals on the right-hand side of (9) and (13) are approximated using the integrands at the grid nodes, considering the areas of the triangles.
Approximating the integrals in (9) and (13) using their discrete analogues leads to a system of linear algebraic equations:
Here, is the total number of triangulation nodes at the external and internal borders and and the factor depends on what type of boundary the node belongs to. If the node is at , then the coefficient is 1; otherwise, it is equal to the value of . The values of represent the coefficient of mutual influence of the points and and are formed by expressing the integrals in a discrete form. Equation (14) can be solved using either direct methods or iterative methods. Both of these methods have been verified to ensure that the solution does not depend on the selected solution method.
The calculation algorithm has been tested using two methods: checking the implementation of the “reciprocity principle” for the four-electrode array and comparing the calculation results with known solutions for a two-layered medium.
4. Reciprocity Principle
Modern ERT equipment can take high density measurements while simultaneously recording data from many measuring electrodes [13, 14]. As a result, more data can be obtained in one experiment. In addition, the same electrodes can act as either measuring or source electrodes in different experiments. Thus, it is possible to diversify the types of measuring electrode arrays without large expenses. To optimize the equipment settings and to accelerate the measurements, the principle of reciprocity is applied in geophysics. Based on this principle, the potential difference between the points , as measured by the installation , does not change if and exchange roles, that is, changing the source electrodes to measuring electrodes and vice versa. Therefore, the principle of reciprocity can be used to reduce the number of measurements by changing the roles of the electrodes.
To test our algorithm, the reciprocity principle has been verified numerically. Let be the numbers of the equidistant positions of electrodes on a measuring line of length . Assume that the source and sink electrodes and occupy positions and , respectively. Then, the numerical experiments for and obtain the values of the potential difference between the electrodes and , , . If electrodes and are located at positions and , , in a series of experiments, then the potential difference between the electrodes with numbers and should be equal to obtained in the previous experiment.
This test was performed for a four-electrode array, exchanging the roles of the source electrodes with those of the measuring electrodes and then calculating the potential difference between points and . The value of was set to . Calculations were performed for two options of the source electrodes and : (a) placed at positions and and (b) placed at positions and . In the first case, the successive positions of the measuring electrodes and were considered. Then, the roles of the electrodes were changed, placing and at points and and then calculating the potential difference between positions and . The average relative difference between the potential differences and did not exceed 2%. Figure 3(a) shows the sequence of , , for the source electrodes located at and . Figure 3(b) plots the results for exchanging the roles of the source electrodes with those of the measuring electrodes
(a)
(b)
The same procedure was executed for the second test for the source electrodes positioned at and . In contrast to the first case, the potential difference between the electrodes and , , was computed excluding the values . The experiments were then performed by exchanging the roles of and .
The relative deviation of the potential difference in the second test did not exceed 2%. Therefore, both tests demonstrate the implementation of the “principle of reciprocity.” The test results of case (b) are shown in Figure 4.
(a)
(b)
Thus, these tests demonstrate the implementation of the principle of reciprocity in our numerical models.
5. Numerical Results and Discussion
A numerical solution of (14) was tested on a known solution for a two-layered medium, which is given, for instance, in [18], and the results are presented in Figures 5 and 6. The calculations were performed for a three-electrode system and a two-layered medium model with a flat surface, a high-resistivity base , and a conductive base .
The influence of the depth of the second layer was confirmed for both decreasing resistivity models and increasing resistivity models. The maximum relative error in calculating the apparent resistivity was not greater than 2%.
This test allows us to determine the admissible dimensions of the computational domain and the mesh thickening necessary to achieve the desired accuracy of approximately 1-2% in the apparent resistivity function.
When the conductivity of the second medium was less than that of the first, the computational region could be smaller than when the second medium had a substantially (~100) larger conductivity. In the calculations, a certain length scale was used as the unit of length, and all the geometric parameters of the model (height of the relief, depth of the second layer, and length of the measuring profile) are expressed in this unit. The symmetry with respect to the abscissa axis is taken into account to halve the number of nodes in the mesh. The resistivity of the first medium is used as the unit of resistivity as well as the units of electrical conductivity. For instance, in the case of , the dimensions of the triangulated oval are equal to 4 by 2. Here, the length of the line along which the mesh gridding is conducted is equal to 2, and the number of nodes on this line varied from 200 to 300. The mesh expanded when removed from the central line. If the expansion coefficient was 8, that is, the size of the largest cell was approximately 8 times greater than the size of the smallest cell, then the parameters are consistent with approximately 16600 triangles. For media with , the dimensions of the computational domain were 9 by 6, with 320 points on the condensation line and a condensation coefficient of 25. The total number of nodes in this case was approximately 10400. The calculation time for the iterative method did not exceed one minute on a personal laptop, whereas a direct method for solving a system of linear equations takes approximately 40 minutes.
Figures 7 and 8 show the solutions of the system of integral equations for the shape of the relief, which is given analytically by the following formula:
In this case, the relief has a shaft shape, as shown in Figure 1, with the underlying layer at a depth of and a conductivity of . The discontinuity at one point of the solution in Figure 7 corresponds to the point of application of the pointwise source electrode. The apparent resistivity, calculated along the profile with a length of 2 in the transverse direction of the relief, is shown in Figure 8.
To determine the effect of the relief on the sounding curves, we consider the results of the numerical simulation for a two-layered medium with a relief of the ground surface. Numerical simulations were performed for reliefs in the form of a shaft (Figure 1) and a valley (Figure 9) with an underlying layer at a depth of . To determine the effect of the resistivity of the underlying layer on the apparent resistivity anomaly associated with the relief, we perform the calculations for varying resistivity values of the underlying layer for models with both conductive and high-resistivity bases.
Figure 10 shows the simulation results for a medium with a conductive underlying layer with the relief of the ground surface in the form of the valley (Figure 9) and without relief (horizontal plane) for varying resistivity values Ωm of the underlying layer. The resistivity of the first layer is set equal to Ωm.
Figure 11 shows the simulation results for the same medium with a higher resistivity base. In this case, the resistivity of the upper layer is Ωm, while that of the underlying layer varies among Ωm.
Figures 12 and 13 show similar results to Figures 10 and 11, but the relief is in the form of a shaft (Figure 1).
The depth of the second layer also influences the sounding curves. To study the influence of the depth of occurrence of the second layer, we perform numerical simulations on a two-layered medium with a relief for different values of . Figure 14 shows the simulation results for m for a decreasing resistivity model, where and Ωm.
For all the cases, the apparent resistivity curves go to the asymptotic values of the second layer . Evident anomalies are shown near the relief irregularities. The vertices of convexity or concavity of the relief form are marked at the sounding curves with a minimum or a maximum, respectively. In the case of a valley relief, the anomalies were expressed more than for a hill relief. In addition, when performing numerical calculations to reach the asymptotic values of the second layer, the computational domain for the model with a high-resistivity base was twice as large as that for the model with a low-resistivity base.
The calculations show that the apparent resistivity anomalies associated with the relief are increasingly expressed when the resistivity of the underlying layer is smaller than that of the upper layer and the lower layer is closer. Conversely, if the underlying layer has a resistivity that is much higher than that of the upper layer, the influence of the relief is less pronounced.
To define the influence of the relief on the inversion results, we calculated the apparent resistivity curves for a relief form similar to that depicted in Figure 1. Numerical simulations were performed for the following parameters: length of measuring line of 235 m; 48 electrodes; depth equal to the height of the relief, namely, 20.9 m; elevation angle 40 degrees; resistivity Ohm·m; resistivity Ohm·m; and = 5 m. The above problems were solved for the three-electrode arrays and . The successive positions of the electrode arrays correspond to a common practice for ERT: the source moves from left to right, and moves in the opposite direction in the steps of MN. The pseudosection obtained for these parameters is shown in Figure 15.
Then, the synthetic data shown in Figure 15 are entered into the 2D resistivity imaging programs Res2DInv [20] and ZondRes2D [21]. The inversion results obtained using Res2DInv and ZondRes2D are depicted in Figures 16 and 17, respectively. As shown in Figures 16 and 17, the relief causes distortions in the interpretation data, especially when using the program Res2DInv. When using the inversion program ZondRes2D and for the considered example, the distortions are less expressed (Figure 17). In fact, the similarity of the original model and the inversion results shown in Figure 17 further verifies our numerical method.
6. Conclusion
Although we discussed the problem of a two-dimensional structure, the electrical field generated in a medium is three-dimensional. The integral equations method has been shown to be highly efficient for this type of field calculation. The calculations show that the apparent resistivity anomalies associated with the relief lead to distortions of the interpretation results in 2D inversion programs. Modelling based on the integral equations method can estimate the level of false anomalies that appear due to topography. Future research will calculate the effect of the relief for more complicated media structure based on field experiments and models. The ultimate goal of this research is to eliminate false anomalies that appear due to topography from the interpretation results.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors express their deep gratitude to Pr. I. Modin and A. Skobelev for allowing them to compare their results from the inversion programs Res2DInv and ZondRes2D. This opportunity was provided during the corresponding author’s visit to Moscow State University. This work is supported by a grant titled “The Best Teacher-2016” and Grant no. 1264/GF4 of the Ministry of Education and Science of the Republic of Kazakhstan.