# Modelling the Influence of Ground Surface Relief on Electric Sounding Curves Using the Integral Equations Method

Authors: Balgaisha Mukanova; Tolkyn Mirgalikyzy; Dilyara Rakisheva
Journal: Mathematical Problems in Engineering
 (2017)

---

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


---

## Body

## 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(A,B) 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σ1 and σ2 (corresponding to specific resistivities ρ1 and ρ2, resp.). Assume that the lower layer is disposed horizontally at depth h 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.Figure 1
Model of a two-layered medium with ground surface relief.In practice, the voltage is applied to the medium through two source electrodes, designatedA and B. However, based on the principle of superposition of electric fields, we consider only the potential created by the point source A.LetГ1 be the ground surface and let Г2 be the contact boundary between the media with conductivities σ1 and σ2. The potential field is the sum of the simple layer potential defined on the boundaries and the potential of point source A in the top boundary Г1:(1)UP=I2πσ1AP+I4πσ1∫Γ1q1MPMdΓM+∫Γ2q2MPMdΓM.The factorI/4πσ1 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 q1 and q2. According to the properties of the simple layer potential, U(P) satisfies Laplace’s equation in the domain Ω1∪Ω2 except at the borders Γ1 and Γ2 and at point A (Figure 1). The function U(P) must satisfy the boundary condition at the inner boundary between two materials Г2 and on the surface of the medium Г1. The first of these conditions enforces the conservation of the current through the contact boundary and is written as(2)σ1∂U∂n1=σ2∂U∂n2.Indicesi=1,2 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:(3)uiP=∫ΓiqiMdΓMPM,i=1,2.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Ω2 is equal to the density of a simple layer multiplied by 4π [19]. If the normal to the boundary Г2 is directed from medium 2 to medium 1, then(4)∂u2∂n2-∂u2∂n1=4πq2P.In this case, the following equalities hold:(5)∂u2∂n2=2πq2P+∫Γ2q2Mcos⁡ψPMPMdΓ2M∂u2∂n1=-2πq2P+∫Γ2q2Mcos⁡ψPMPMdΓ2M.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Γ2is written as(6)σ2∂u2∂n2=σ1∂u2∂n1-σ2-σ14πσ1I∂U0∂n+∂u1∂nP.Substituting (5) into (6), we obtain(7)2πq2P=σ1-σ2σ2+σ1∫Γ2q2Mcos⁡ψPMPMdΓ2M+4πσ1I∂U0∂n+∂u1∂nP.Hence, the following integral equation is obtained:(8)q2P=λ2π∫Γ2q2Mcos⁡ψPMPMdΓ2M+λ2π4πσ1I∂U0∂n+∂u1∂nP.Here,λ=(σ1-σ2)/(σ1+σ2) denotes a reflection coefficient at the boundary between media 1 and 2. Substituting the potential for source electrode A and u1(P) into (8), we obtain (9)q2P=λ2π∫Γ2q2Mcos⁡ψPMPMdΓ2M+λ2π∫Γ1q1Mcos⁡ψPMPMdΓ1M+λπ∂∂n1AP,P∈Γ2.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:(10)I4π∂u1∂nPΓ1+∂u2∂nPΓ1+σ1∂U0∂nΓ1=I2πδr→.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 Г1. 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(11)I4π∂u1∂nPΓ1+∂u2∂nPΓ1+σ1∂U0∂nΓ1=I2πδr→.Substituting (11) into (10), we obtain(12)2πq1P+∫Γ1q1Mcos⁡ψPMPMdΓ1M+∂u2∂nPΓ1+4πσ1I∂U0∂nΓ1=2δr→.The above equation yields another integral equation:(13)q1P=-12π∫Γ1q1Mcos⁡ψPMPMdΓ1M-12π∫Γ2q2Mcos⁡ψPMPMdΓ2M+1πδr→-∂∂n1APΓ1,P∈Γ1.Thus, we have two integral relations, (9) and (13), for the two unknown functions of the simple layer densities q1(P) and q2(P), which are defined at the ground surface Г1 and at the contact boundary Г2, 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 Г1 and Г2 by their finite parts representing the oval-shaped regions. The boundary Г1 is mapped on the plane onto which the irregular triangular mesh is constructed. The grid is condensed along the Ox 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 Figure2, and the unknown functions q1(P) and q2(P) 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.Figure 2
Typical triangulation of the computed region (the positivez-direction is directed downwards).Approximating the integrals in (9) and (13) using their discrete analogues leads to a system of linear algebraic equations:(14)qi=ki2π∑j=1NqjAi,j+kiπBi,i=1,N-.Here,N is the total number of triangulation nodes at the external and internal borders Г1 and Г2 and the factor ki depends on what type of boundary the node belongs to. If the node is at Г1, then the coefficient is 1; otherwise, it is equal to the value of λ. The values of Aij represent the coefficient of mutual influence of the points i and j 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 arrayAMNB 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 MN, as measured by the installation AMNB, does not change if AB and MN 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. Letk=0,…,K be the numbers of the equidistant positions of electrodes on a measuring line of length L. Assume that the source and sink electrodes A and B occupy positions k=0 and k=1, respectively. Then, the numerical experiments for A and B obtain the values of the potential difference ΔUij between the electrodes Mi and Mj, i=2,…,K-1, j=2,…,K-1. If electrodes A and B are located at positions Mi and Mj, i=2,…,K-1, j=2,…,K-1 in a series of experiments, then the potential difference between the electrodes with numbers 0 and 1 should be equal to ΔUij obtained in the previous experiment.This test was performed for a four-electrode array, exchanging the roles of the source electrodes(A,B) with those of the measuring electrodes (M,N) and then calculating the potential difference between points M and N. The value of K was set to K=25. Calculations were performed for two options of the source electrodes A and B: (a) placed at positions 0 and 1 and (b) placed at positions 4 and 5. In the first case, the successive positions of the measuring electrodes Mi and Mi+1 were considered. Then, the roles of the electrodes were changed, placing A and B at points Mi and Mi+1 and then calculating the potential difference ΔU01 between positions 0 and 1. The average relative difference between the potential differences ΔU01 and ΔUi,i+1 did not exceed 2%. Figure 3(a) shows the sequence of ΔUii+1, i=2,…,K-1, for the source electrodes located at 0 and 1. Figure 3(b) plots the results for exchanging the roles of the source electrodes (A,B) with those of the measuring electrodes (M,N)Figure 3
Potential differenceUi-Ui+1, depending on the position of MN (a) and the potential difference U0-U1 after exchanging the roles of AB and MN with A at i and B at i+1 (b).
(a)
(b)The same procedure was executed for the second test for the source electrodes positioned at4 and 5. In contrast to the first case, the potential difference between the electrodes M0 and Mi, i=1,…,K, was computed excluding the values i=4,5. The experiments were then performed by exchanging the roles of AB and MN.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 Figure4.Figure 4
Difference in potentialU0-Ui depending on the positions of N (a) and the difference of potential U4–U5 after exchanging the roles of AB and MN, with A placed at 0 and B placed at i (b).
(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 AMN and a two-layered medium model with a flat surface, a high-resistivity base (ρ1<ρ2), and a conductive base (ρ1>ρ2).Figure 5
The function of the apparent resistivity for a model with a conductive base: (1) a solution built in [18], (2) a solution obtained by solving equation (14).Figure 6
The function of the apparent resistivity for a model with a high-resistivity base: (1) a solution built in [18]; (2) a solution obtained by solving (14).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σ2=10σ1, 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 σ2=0.001σ1, 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.Figures7 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:(15)zx=0.06cos⁡2πx+0.25-1,-0.25≤x≤0.750,x∉-0.25,0.75.Figure 7
Distribution of the density of a simple layerq1(P) at the surface Г1.Figure 8
Distribution of the densityq2(P) at the surface Г2.In this case, the relief has a shaft shape, as shown in Figure1, with the underlying layer at a depth of z=0.5 and a conductivity of σ2=10σ1. 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 (Figure1) and a valley (Figure 9) with an underlying layer at a depth of h=0.5. To determine the effect of the resistivity of the underlying layer on the apparent resistivity (ρa) 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 9
Model of a two-layered medium with relief in the form of a valley.Figure10 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 ρ2=10,30,50 Ωm of the underlying layer. The resistivity of the first layer is set equal to ρ1=100 Ωm.Figure 10
The shape of the surface relief and the apparent resistivity curves: (1) for a two-layered medium model with a relief; (2) for a model of a two-layered medium with a horizontally plane surface.Figure11 shows the simulation results for the same medium with a higher resistivity base. In this case, the resistivity of the upper layer is ρ1=10 Ωm, while that of the underlying layer varies among ρ2=30,50,70 Ωm.Figure 11
Surface relief shape and apparent resistivity curves: (1) for a two-layered medium model with a relief; (2) for a model of a two-layered medium with a horizontally flat surface.Figures12 and 13 show similar results to Figures 10 and 11, but the relief is in the form of a shaft (Figure 1).Figure 12
Surface relief shape and apparent resistivity curves: (1) for a two-layered medium with a horizontally flat surface; (2) for a two-layered medium with a relief.Figure 13
The shape of the surface relief and the apparent resistivity curves: (1) for a two-layered medium with a relief; (2) for a model of a two-layered medium with a horizontally plane surface.The depth of the second layerh 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 h. Figure 14 shows the simulation results for h=0.25,0.5,0.75 m for a decreasing resistivity model, where ρ1=100 and ρ2=10Ωm.Figure 14
The influence of the occurrence depth of the second layer on the apparent resistivity curves.For all the cases, the apparent resistivity curves go to the asymptotic values of the second layerρ2. 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 Figure1. Numerical simulations were performed for the following parameters: length of measuring line AB of 235 m; 48 electrodes; depth h equal to the height of the relief, namely, 20.9 m; elevation angle 40 degrees; resistivity ρ1=100 Ohm·m; resistivity ρ2=10 Ohm·m; and MN = 5 m. The above problems were solved for the three-electrode arrays AMN and BMN. The successive positions of the electrode arrays correspond to a common practice for ERT: the source A moves from left to right, and B moves in the opposite direction in the steps of MN. The pseudosection obtained for these parameters is shown in Figure 15.Figure 15
Pseudosection calculated for a model of a two-layered medium with ground surface relief.Then, the synthetic data shown in Figure15 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.Figure 16
Interpretation data obtained by the program Res2DInv for the model of a two-layered medium with ground surface relief.Figure 17
Interpretation data obtained by the program ZondRes2D for the model of a two-layered medium with ground surface relief.


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



---

*Source: 9079475-2017-07-26.xml*