Abstract
The direct task of the subsurface exploration of a homogeneous medium with surface relief by the resistivity method is analyzed. To calculate the resistivity field for such a medium, the method of integral equations was successfully applied for the first time. The corresponding integral equation for the density of secondary current sources on the surface of the medium was established. The method of computational grid construction, adapted to the characteristics of the surface relief, was developed for the numerical solution of the integral equation. This method enables the calculation of the resistivity field of a point source on a surface that is not smooth and allows for steep ledges. Numerical examples of the calculation of resistivity fields and apparent resistivity are shown. The anomalies of apparent resistivity arising from the deviation of the surface shape from a flat medium were quantitatively established as model examples. Calculations of apparent resistivity for the direct current sounding method were carried out using modifications of the electrical tomography approach.
1. Introduction
Accounting for the impact of ground surface relief during the solution of direct electrical problems is an important issue in the field of electrical resistivity tomography. Currently, the most commonly used methods for solving direct and inverse problems are grid methods (e.g., finite difference method and finite element method) [1–6]. Within these grid methods, several approaches to account for the influence of surface relief were developed.
In the finite difference method, two basic algorithms accounting for the effects of surface topography were developed. The first and most widely used algorithm is the geometric transformation of cells (node position) in the grid in accordance with the relief of the ground surface [7]. The second method of accounting for ground surface topography for the finite difference method is the approximation of the transition from a five-point cross pattern to a triangular pattern, which allows the use of triangular cells, similar to cells used in the finite element method [8, 9]. It is more efficient to transition to triangular cells (generalized finite difference method) in combination with a specific approach to the elliptic problem, in which the task is carried out along the borders of the grid lines [10]. Clearly, all the circuits of the cell that have a rectangular shape contribute to geometric error when considering the effect of surface topography. It is possible to reduce the error only by increasing the number of cells, which leads to a significant increase in computational cost.
The finite element method, in conjunction with the direct problem solution methods of direct current, is rapidly developing. This applies to both the configuration of the grid and the mathematical apparatus. In the finite element method, there are several approaches to accounting for surface topography in solving direct problems [11–13]. As with the finite difference method, the first method of accounting for surface topography is a transformation of cells. The second approach is to construct a special irregular grid depending on the form of the topography [14] and the degree of heterogeneity of the section. It should be noted that BERT systems [15] use these approaches.
Finite difference methods for modeling direct current problems are based on the numerical solution of an elliptical equation using finite difference schemes:
is the extraneous current flowing into the volume unit of land , is a conductivity of the medium, and is a potential.
The idea of the method is to divide the modeled field into a grid of rectangular blocks, each of which is characterized by its value of resistivity. In the boundary nodes of the grid, the solution of (1) must satisfy the boundary conditions, including Dirichlet (the field must be equal to zero at large distances from the source, in the lower and side corners of the grid), Neumann (the normal component of the field is equal to zero at the ground surface boundary), or mixed boundary conditions:where is the angle between (the radius-vector from the source to the node) and (outward normal to the boundary).
The idea of the finite element method consists of dividing cells into triangular shapes (elements) and presenting potential in the form ofwhere is a set of basic functions. In this form, the desired function for each of the elements, based on the values of its nodes, is approximated. Equations for minimizing misties are used for each node and, together with the equations for all other elements of the grid, are used to form a general matrix equation. Solving this equation, the desired potential is defined. The boundary conditions for the boundary elements of the grid are the same as in the finite difference method.
They are universal and useful for many practical purposes but have several disadvantages. In particular, these disadvantages are as follows:(1)low accuracy of the field’s first derivative computation and low accuracy of subsequent computations of apparent resistivity;(2)problems with producing artificial boundary conditions on boundaries of the computation domain and big distortions of the apparent resistivity [16];(3)problems with creating the grid, resulting from the complicated form of the boundary of the computation domain.
On the other hand, purposes of geophysics have high requirements on accuracy and speed of field computations. Therefore, to control the operation of grid methods, it is appropriate to apply other methods based on different ideas that enable the solution of the problem related to the direct method of mapping electrical resistivity.
As an alternative method of solving the direct problems, satisfying both of these requirements is the method of integral equations. The main new feature of the presented work is that we apply the method of integral equations for calculation of fields in a medium with a relief surface for the first time. Different applications of the integral equations method in geophysics are considered in articles [17–24] and are well established among the methods of mapping electrical resistivity.
Advantages of the method of integral equations in comparison with grid methods are the following:(i)A three-dimensional problem is reduced to two-dimensional one, which significantly decreases the calculations.(ii)It has a higher accuracy of calculations of an electrical field and apparent resistivity.(iii)It allows taking into account the shape of the relief and the contact boundaries of the medium almost without distortions.(iv)It has a clear physical meaning.
The main disadvantage of the method is that it is intended for a special class of models of mediums. Namely, it is required that medium has sharply defined geoelectrical boundaries, and it is not as versatile as the grid methods. However, because of its specialization, the method is highly efficient for its class of mediums’ models.
2. The Problem with the Effect of Ground Surface Relief
The electrical exploration method is a method of excitation and observation of the electric field in order to obtain information about the structure of a geological section. One of the main methods of modern electrical exploration is vertical electric sounding with the modification of 2D electrical tomography.
Essentially, electrical tomography uses a chain of grounded metal electrodes on the ground surface, placed along the observation profile at equal distances from each other. At a certain point, an electrical voltage is applied on two electrodes, which are included in the mentioned system. Connections of certain pairs of electrodes are held alternately, according to the protocol. Voltage on the source electrodes excites an electric current in a conducting inhomogeneous medium. In turn, this current produces electric fields in the subsurface medium, and the spatial and amplitude characteristics are strongly dependent on factors such as the structure of the geologic section, subsurface heterogeneity, and ground surface relief [8–11, 13, 14, 25].
Ground surface relief greatly complicates the task of solving the inverse problem of finding the distribution of according to physical measurement of the components of the electric field (, ) on the ground surface. In almost all methods of electrical prospecting, the measurements are influenced by the ground surface topography due to a near-surface heterogeneity, which causes a change in the flow lines. In [11] the general effects of topographic relief for a remote current source are discussed. In this case, the current flow lines diverge on the convex part of the surface and converge on the concave part. As a result, the directions of the electric field lines change. Because of this effect, the concave shape increases and the convex shape of the surface decreases the value of the apparent resistivity.
The complexity in the description of the resistivity field near the corner points of the ground surface should be noted. The disadvantages of these effects are more important when we have sharp boundaries and sharp corners in the model because even their insignificant smoothing will significantly change the field of apparent resistivity.
The overlay of distorting effects caused by the inhomogeneity of ground surface topography creates a complex picture of the resulting geoelectric pseudosections. Excluding these distorting factors in the solution of the inverse problem, or some version of a two-dimensional or three-dimensional inversion, it is impossible to get closer to the construction of final sections that are reasonably accurate.
3. Application of the Method of Integral Equations to the Solution of the Direct Problem of Electrical Prospecting with Direct Current under Conditions of Ground Surface Relief
3.1. The Method of Integral Equations and Its Development Based on Ground Surface Relief
The idea of this method is to provide the electric field as a sum of the primary field (generated by current electrodes, in this case, by one electrode) and the field of secondary charges (occurring when the electric current flows in points of homogeneity violation of the subsurface medium and on the surface of the medium):where is the vector of the primary electric field and is the vector of total electric field of the secondary charges.
Contact boundaries and heterogeneous inclusions of the geoelectric section act as secondary creators of the electric field. The field computation problem is reduced to the system of integral equations on the current density of secondary sources, inducted on contact surfaces of conductive mediums and on the ground surface of the medium. The mathematical description of this event leads to the Fredholm equations of type II, which, after discretization of the problem, leads to a system of linear algebraic equations (SLAE).
Let us consider a three-array of Schlumberger , with one source () and two receiver () electrodes. Source electrode is located in the center of the coordinates (Figure 1).
The electrostatic potential at points in the medium satisfies the differential equation of Laplace for stationary fields. The condition of decreasing the potential at infinity and the condition at the ground surface boundary must be satisfied:
Compared with a case when the surface is plain, the presence of ground surface relief makes a problem formulation for the method of integral equations more complicated. This is because we can use the method of reflections at a plain surface of discontinuity at the ground surface, which allows us to use a relatively simple system of integral equations [23, 24] (in the case of multiple internal contacts of boundaries) [22]. To calculate the field at the ground surface, we define a distribution of current densities at the surface, which compensates for an external source to satisfy the physical boundary condition without current normal to the ground surface boundary:
Here, is current in a direction normal to the surface, which might arise from a source electrode in the absence of compensatory secondary currents. According to Ohm’s law (6), the condition can also be written in terms of the field strength of the external source:
Here, is the normal to the surface of the component of the electric field created by the source electrode.
We write the integral equation for the distribution of secondary sources using boundary condition (7). Let be a point on the surface of the medium , denoted by , the current density of the secondary source at point (Green’s function of a point source); the mutual influence of all other sources on the source of is calculated under the integrand. From the integral equation, it is required to calculate the current density of the secondary source, :where is the normal of the component of the electric field from the primary source at the point ; is the potential of a point source placed at the surface point generated at point (the source point of the Green function); is the electric current density of the secondary source at point ; and is the area of the surface element at point . Thus, integral equation (8) expresses the physical fact that the points of the total surface current are equal to zero, except for the point of application of the source electrode.
3.2. The Transition from the Integral Equation to the System of Linear Algebraic Equations
The solution of integral equation (8) can be obtained by transforming the system of linear algebraic equations. To describe the distribution of current density on the surface, we will apply a grid with a logarithmic scale expanding radially, adapted for the source electrode position and the shape of the surface relief (Appendix). Thus, the surface is divided into elementary cells with vertices at grid nodes.
Discrete approximation of integral equation (8) can be written in system form aswhere is the number of cells broken on the surface ; are the number of cells; is the cell with index ; is the required value of the current density of secondary sources; and is the value of the electric field from the primary sources. determines the influence of coefficients at the th points on the power in the th point. is calculated by numerical integration at the normal to the surface at of the derivative of Green’s function in the coordinates of the sources, located within cell . Point is the geometric center of cell . Function (10) can be written in the following way:where is the geometric center of cell and is area of cell .
In our case (homogeneous half-space with conductivity ), , it is calculated aswhere is the vector connecting points and .
Thus, the coefficients of linear algebraic equations (9) are calculated and the system can be written in matrix formwhere is a square matrix of coefficients of the system, is a vector of unknown values of the density of secondary sources, and is a vector of independent terms (right parts) of (9).
3.3. Stages of Solving a Problem by the Method of Integral Equation
(1) Creating the Model and Grid Construction. The medium is described, and the function of the surface relief is determined. The grid is constructed to break the boundaries of the cells approximating the secondary current source and the parameters of each cell (number, location of the geometric center, and area) are defined.
Construction of a grid and determination of the geometry of the cells are described in detail in the Appendix. To calculate the total current flowing through the surface of the cells, it is necessary to calculate the area for each cell. The area of the cell is identified by the following formula:where , are low and high angles and , are top and bottom radii of the cell with the number .
(2) Construction of the Matrix with Coefficients of Mutual Influence. The coefficient matrix, the vector of unknown values, and vector of independent terms of system (13) are identified as follows:where here,
, , is a vector directed from point to point ; , is the value of the current flowing from the current source electrode, and is the vector directed from point , where the source electrode is placed, to the point with the number .
(3) Solving the System of Linear Algebraic Equations and Determining the Density of the Secondary Current Sources . System (15) is fulfilled and can be solved by the elimination method of Gauss or any iterative methods.
(4) Calculating the Potential on the Surface. After solving the system of equations and finding the distribution of secondary sources, the values of the potentials on the surface are calculated by use of numerical integration. Below, the formula for calculating the potential of the secondary source at point is presented:
Here, is the potential of the source electrode at point .
(5) Calculating the Functions of Apparent Resistivity. Calculating the functions of apparent resistivity along the surface is defined by the relation where is the current source electrode and is the geometric coefficient for the three-electrode Schlumberger array.
4. Numerical Results and Discussion
For a full explanation of the anomalies due to surface relief, the electric current density of the surface (with secondary sources for different surface shapes) computed by the integral equation (8) is shown in Figures 2(a), 2(b), and 2(c).
(a)
(b)
(c)
The results obtained from the numerical simulations are presented in Figures 3–7. Calculations were made for a three-electrode Schlumberger array for several models of the Earth’s surface relief for a medium with resistivity Ωm.
(a)
(b)
(a)
(b)
Here, we examine the simulation results for the model with a positive form of ground surface relief in the form of a hemispherical (semicircular) crown with sharp angles (Figure 3) and also for the same model results with smooth angles.
From these curves of apparent resistivity, the anomalies can be observed in the vicinity of the surface irregularities, especially at the corners. On the flat areas of the surface, the of the curves are close to the resistivity of the medium. As you approach the corner point, there is a rise and a sharp increase in , and the uppermost minimum apparent resistivity is at the vertex. Then, is found to increase in the vicinity of the next corner point. This is followed by a smooth decrease in as it approaches the resistivity of the medium. From the results of the simulation of the curves of for ground surface forms with smooth angles of inclination, it can be seen that the magnitude of the anomaly decreases with smoothing of the corner points. Increased races are explained by distancing receiver electrodes from the current source electrode.
Figure 4 presents the results of the simulation model with the surface in the form of a hemispherical (semicircular) recess. In the corner points of the surface, a decrease in is observed, and the center of the recess corresponds to the maximum. For the surface model with smooth angles, the anomaly is less apparent. It should be noted that, at the corners of the hemispherical (semicircular) recess surface, the anomaly is significantly smaller than for the hemispherical (semicircular) convex surface.
Here, we examine the results of the simulation model in which there are two shapes of the ground surface: a positive (hemispherical convexity) and a negative (hemispherical recess) relief (Figure 5). For this model, there are similar patterns of behavior of the curves, as with previous models. The value of increases near the top of the concave corner of the curve and decreases near the vertices of the convex corner. Additionally, there is a smooth transition from a minimum to a maximum corresponding to the transition from the positive to the negative forms of the surface.
The simulated values of models where the surface is in the form of a symmetrical triangle, with positive and negative forms, are shown in Figures 6 and 7. The calculations are made for two slope angles: , 22.5°. For positive relief, increases in the angles of inclination result in a decrease in found in the extreme corner points of the curve. For negative relief, the opposite behavior is observed. It should be noted that the magnitude of the anomalies increases with increasing angle of slope for either positive or negative reliefs in the symmetrical triangles. For the same angles of slopes and for positive relief, anomalies are greater than those for negative relief.
5. Conclusion
This paper presents a mathematical model that allowed us to construct a function of apparent resistivity that is constructed using the basic laws of formation of the electric field associated with the influence of ground surface relief. The numerical mathematical algorithm and computer program solve the direct problem of electrical prospecting with direct current for arbitrary two-dimensional and three-dimensional media and account for the effects of ground surface relief. The results will improve the quality of work and the development of computer technology used in geophysics.
This research can be continued in several different directions. Firstly, the method should be adapted and tested on real formats of surfaces obtained in the field measurements. Secondly, it is necessary to generalize the method to the case of more complex models of medium containing local inclusions, buried relief, or embossed layered structures. The last option is the most common for geophysical research. Further, it is necessary to develop a quantitative methodology of accounting of anomalies associated with the relief, in the interpretation of data. This problem requires the formulation and solution of inverse problem in terms of the integral equations method.
Appendix
Construction of a Computational Grid on the Ground Surface Relief
One of the important steps in solving direct problems of electrical exploration with direct current is to build a computational grid. The quality of the computational grid affects the accuracy of the results, the convergence of the iterative method, and time of the decision of the whole problem.
In the method of integral equations the grid is constructed of the cells broken on the surface, approximating the secondary current source.
To describe the distribution of current density on the surface, we create a grid with a logarithmic scale expanding radially, adapted for a source electrode position and the shape of the surface relief. In [26], the comparison of the results of calculations of the potential for homogenous and inhomogeneous grids with a logarithmic scale expanding is presented. It is shown that, in calculations based on inhomogeneous grids condensing in the vicinity of the source electrode, the relative accuracy of the calculation is higher than on a homogeneous grid. Initially the grid nodes are defined in the polar coordinate system associated with the surface (Figure 8). This method allows one to parameterize the complex structured relief including vertical plumb lines and even recess.
The initial values and steps for each coordinate and are set. Then the coordinates are determined. Step of the corner is defined by a constant , and the step of the radius is determined by the logarithmic expanding. Logarithmic expanding step along the radius is calculated by the change of variables in the domain of integration : . The results obtained in the polar coordinate system are converted into Cartesian and the coordinates of grid points are calculated by taking into account the surface relief . Coordinates are found by solving the differential equation by Runge-Kutta method; the coordinates are determined by the function of relief . In Figure 9 a grid is constructed by the described algorithm for the relief function when and when , .
Incremental algorithm of constructing a grid on the surface of the relief:(1)determination of the relief function ;(2)calculation of the polar grid nodes , , :(3)conversion of coordinates obtained in the polar coordinate system to Cartesian coordinate system :(4)calculation of the coordinate functions of the relief :
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors would like to thank the anonymous referees whose comments aided to improve the paper. This work is supported by the Ministry of Education and Science of the Republic of Kazakhstan, Grant no. 46-12.02.2015.