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

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


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][4][5][6][7][8][9][10][11][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][10][11][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][10][11][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.

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 ℎ 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 D 1 be the ground surface and let D 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  in the top boundary D 1 : The factor /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  1 and  2 .According to the properties of the simple layer potential, () satisfies Laplace's equation in the domain Ω 1 ∪ Ω 2 except at the borders Γ 1 and Γ 2 and at point  (Figure 1).The function () must satisfy the boundary condition at the inner boundary between two materials D 2 and on the surface of the medium D 1 .The first of these conditions enforces the conservation of the current through the contact boundary and is written as Indices  = 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: 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 D 2 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 Γ 2 is written as Substituting ( 5) into (6), we obtain Hence, the following integral equation is obtained: Here,  = ( 1 −  2 )/( 1 +  2 ) denotes a reflection coefficient at the boundary between media 1 and 2. Substituting the potential for source electrode  and  1 () 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 D 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 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  1 () and  2 (), which are defined at the ground surface D 1 and at the contact boundary D 2 , respectively.

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 D 1 and D 2 by their finite parts representing the oval-shaped regions.The boundary D 1 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  1 () and  2 () 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 D node is at D 1 , 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.

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  = 0, . . .,  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  = 0 and  = 1, respectively.Then, the numerical experiments for  and  obtain the values of the potential difference Δ  between the electrodes   and   ,  = 2, . . .,  − 1,  = 2, . . .,  − 1.If electrodes  and  are located at positions   and   ,  = 2, . . ., −1,  = 2, . . ., − 1 in a series of experiments, then the potential difference between the electrodes with numbers 0 and 1 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  = 25.Calculations were performed for two options of the source electrodes  and : (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   and  +1 were considered.Then, the roles of the electrodes were changed, placing  and  at points   and  +1 and then calculating the potential difference Δ 01 between positions 0 and 1.The average relative difference between the potential differences Δ 01 and Δ ,+1 did not exceed 2%. Figure 3(a) shows the sequence of Δ +1 ,  = 2, . . .,  − 1, for the source electrodes located at 0 and 1.The same procedure was executed for the second test for the source electrodes positioned at 4 and 5.In contrast to the first case, the potential difference between the electrodes  0 and   ,  = 1, . . ., , was computed excluding the values  = 4, 5.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.
Thus, these tests demonstrate the implementation of the principle of reciprocity in our numerical models.

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 ( 1 <  2 ), 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 sions 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. 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  = 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 (Figure 1) and a valley (Figure 9) with an underlying layer at a depth of ℎ = 0.5.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  2 = 10, 30, 50 Ωm of the underlying layer.The resistivity of the first layer is set equal to  1 = 100 Ω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  1 = 10 Ωm, while that of the underlying layer varies among  2 = 30, 50, 70 Ω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 ℎ = 0.25, 0.5, 0.75 m for a decreasing resistivity model, where  1 = 100 and  2 = 10 Ωm.
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 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  1 = 100 Ohm⋅m; resistivity  2 = 10 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.

Conclusion
Although we discussed the problem of a two-dimensional structure, the electrical field generated in a medium is threedimensional.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.

Figure 1 :
Figure 1: Model of a two-layered medium with ground surface relief.

Figure 2 :
Figure 2: Typical triangulation of the computed region (the positive -direction is directed downwards).

Figure 3 :
Figure 3: Potential difference   −  +1 , depending on the position of  (a) and the potential difference  0 −  1 after exchanging the roles of  and  with  at  and  at  + 1 (b).

Figure 3 (
b) plots the results for exchanging the roles of the source electrodes (, ) with those of the measuring electrodes (, )

Figure 4 :
Figure 4: Difference in potential  0 −   depending on the positions of  (a) and the difference of potential  4 - 5 after exchanging the roles of  and , with  placed at 0 and  placed at  (b).

Figure 5 :
Figure5: 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 :
Figure6: 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).

Figure 7 :Figure 8 :Figure 9 :
Figure 7: Distribution of the density of a simple layer  1 () at the surface D 1 .

FigureFigure 11 :
FigureThe 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.

FigureFigure 13 :
Figure Surface relief shape and apparent curves: (1) for a two-layered medium with a horizontally flat surface; (2) for a twolayered medium with a relief.

Figure 14 :
Figure 14: The influence of the occurrence depth of second layer on the apparent resistivity curves.

Figure 15 :
Figure 15: Pseudosection calculated for a model of a two-layered medium with ground surface relief.

Figure 16 :
Figure 16: Interpretation data obtained by the program Res2DInv for the model of a two-layered medium with ground surface relief.