## Abstract

A major difficulty of electrical resistivity forward modelling is the singularity of the potential occurring at the source location. To avoid large numerical errors, the potential is split into a primary part containing the singularity and a secondary part. The primary potential is defined analytically for flat topography, and is classically computed numerically in the presence of topography: in that case, an accurate solution requires expensive computations. We propose to define the primary potential as the analytic solution valid for a homogeneous model and flat topography, and to modify accordingly the free surface boundary condition for the secondary potential, such that the overall potential still satisfies the Poisson equation. The modified singularity removal technique thus remains fully efficient for any acquisition geometries, without any additional numerical computation, and also applicable in the presence of a buried cavity. This approach is implemented with the generalized finite difference method developed on unstructured meshes and validated through the comparison with analytical solutions. Finally, we illustrate in simple 2-D and 3-D cases how the potential depends on the shape of the topography and on the electrode positions.

## INTRODUCTION

Direct current resistivity imaging is a widely used geophysical method. Its range of applications includes mineral resource exploration, groundwater, environmental and geotechnical surveys. 2-D and 3-D investigations are now carried out on a regular basis. In such multidimensional cases, interpretation of resistivity data needs an accurate modelling technique. Field measurements are often performed on a non-flat surface. To discriminate the effect of topography from subsurface inhomogeneities, it is necessary to include topography in both the modelling and the inversion processes.

The ability of taking topography into account in a forward solver greatly depends on the numerical method. The finite difference approach used to be a popular method in the past but restricted to flat topography (Mufti 1976; Dey & Morrisson 1979; Spitzer 1995). Alternatively, the integral equation method is usually valid for simple structures (Dieter *et al.*1969; Okabe 1981). Including complex topography and substructures in the modelling remains difficult with such methods. The finite element approach was first introduced in the resistivity problem by Coggon (1971) and later further developed by many other authors (Zhou & Greenhalgh 2001; Rücker *et al.*2006). It appears to be a flexible method to model both topography and possibly complicated geometries. Initially, finite element modelling codes including topography were defined on distorted grids, using quadrangular or hexahedral elements, and were thus limited to undulating topography (Yi *et al.*2001; Zhou & Greenhalgh 2001). These approaches using block-oriented meshes are also limited regarding the possibility of effective mesh refinements. Only recently, codes have been developed on unstructured tetrahedral meshes (Rücker *et al.*2006; Blome *et al.*2009; Ren & Tang 2010), allowing for complex geometries and local mesh refinements, especially around source locations. In parallel, surface integral methods are also available, with the advantage of only discretizing the top interface (Poullikkas *et al.*1998; Tsai *et al.*2006; Lin *et al.*2011).

Independent of the numerical scheme, a resistivity forward solver consists of solving the electrical boundary value problem for every source location. The solution of such a system of linear equations shows a singularity at the source location. This leads to significant numerical errors in the vicinity of the electrode position. The singularity removal technique introduced by Lowry *et al.* (1989) splits the potential into a singular part, related to the source, and a residual part taking into account a possible non-homogeneous conductivity model. The singular part is defined analytically for flat topography. For non-flat topography, the singular component is usually computed numerically (Blome *et al.*2009). Even with local mesh refinement around the sources, a highly accurate solution is however difficult to be obtained (Rücker *et al.*2006) as the objective is to estimate a singular function. Note that it is carried out only once in an inversion process (Günther *et al.*2006).

We propose a modification of the singularity removal technique to better handle the topography. The basic idea is to use the analytical solution defined for flat topography and homogeneous models as the primary potential and to modify the free surface conditions for the computation of the secondary potential. We discuss the value of the normal derivative of the potential and propose a solution to limit its maximum value. We implement this approach using the generalized finite difference method (GFDM). This method is an extension of the classical finite difference approach on unstructured meshes. Several 2-D and 3-D modelling examples are presented to validate this approach and to illustrate the importance of accurately considering topography.

## SINGULARITY REMOVAL TECHNIQUE WITH TOPOGRAPHY

We first recall the basic resistivity equations. Let us consider a Cartesian coordinate system **x** = (*x*, *y*, *z*). A current source *I* is located at a source point **x _{s}** = (

*x*,

_{s}*y*,

_{s}*z*) in a 3-D conductivity model σ(

_{s}**x**). The electrical potential

*V*(

**x**) is governed by the Poisson equation and the following boundary conditions

_{n}

*V*at the free surface Γ means that equipotentials are perpendicular to the top interface, or equivalently that the current lines are parallel to this surface. Dey & Morrisson (1979) proposed approximate boundary conditions for the external boundaries ∂Ω. Considering that the potential decreases as 1/

*r*, where

*r*is the distance between the source

**x**and a point in the subsurface

_{s}**x**, they proposed the following mixed condition

For a 3-D model invariant along the *y*-direction and a point source at (*x*_{s}, 0, *z*_{s}), the 3-D solution can be obtained with a series of 2-D numerical simulations. A Fourier cosine transform is applied to the potential (Coggon 1971). The modified potential is

*k*

_{y}is the wavenumber related to

*y*. The transformed potential is governed by the Helmholtz equation

*k*

_{y}parameters. The 3-D potential is then obtained by performing the inverse Fourier transform

*k*

_{y}values for a correct modelling are discussed in Queralt

*et al.*(1991).

### Flat topography case

The electrical boundary value problem (eqs 1–3) shows a singularity at the source location. This leads to large numerical errors if *V* is solved directly. A more efficient approach introduced by Coggon (1971) and later improved by Lowry *et al.* (1989) consists of splitting the potential into a primary potential *V*_{p} and a secondary potential *V*_{s}. The primary potential *V*_{p} contains the singularity due to the current source, whereas the secondary potential *V*_{s} contains the effects of a heterogeneous conductivity model. By definition, the total potential *V* is defined as *V* = *V*_{p} + *V*_{s}. Eqs (1) to (2) are rewritten as

*V*=

_{p}*I*/2πσ

_{0}‖

**x**−

**x**‖ is the analytic solution valid for a homogeneous conductivity model and a source located at the surface. For mixed boundary conditions, eq. (10) becomes ∂

_{s}_{n}

*V*

_{s}+ cos θ/

*rV*

_{s}= −∂

_{n}

*V*

_{p}− cos θ/

*rV*

_{p}. The constant σ

_{0}is defined as the conductivity at the source location. In this case, the boundary conditions at the surface Γ are the same for the primary and secondary potential, as the primary potential satisfies ∂

_{n}

*V*

_{p}= 0 at the ground surface. Overall, the Neumann conditions for the total potential

*V*(eq. 2) still holds.

### Arbitrary topography case

In the case of an arbitrary topography, the classical approach consists of solving eqs (8) and (9) and to determine the primary potential satisfying ∂*V*_{p} = 0 (Lowry *et al.*1989). We propose to modify the boundary conditions (eq. 9). For an arbitrary topography case, the total potential is still split into two parts. The primary potential *V*_{p} = *I*/2πσ_{0}*r* remains identical. The only modification resides in the new definition of the Neumann conditions for the secondary potential

**n**is the outward normal at the surface. As above, for mixed boundary conditions, eq. (13) becomes ∂

_{n}

*V*

_{s}+ cos θ/

*rV*

_{s}= −∂

_{n}

*V*

_{p}− cos θ/

*rV*

_{p}. By definition, ∂

_{n}

*V*

_{s}=−∂

_{n}

*V*

_{p}ensures that ∂

_{n}

*V*= 0.

For flat topography, ∂_{n}*V*_{p} = 0 and thus eqs (9) and (12) are identical. For an arbitrary topography, *V*_{p} still contains the singular part but is analytically known. *V*_{s} is function of the topography and of the possibly heterogeneous conductivity model. We extend the approach proposed by Lowry *et al.* (1989) by simply modifying the Neumann conditions for the secondary potential. We now discuss the shape of the secondary potential.

### Singularity in the secondary potential

The singularity in the secondary potential may come from a singular source term (eq. 11) or from singular Neumann conditions (eq. 12). As explained in Appendix A, the source term is not singular if the conductivity is locally homogeneous around the electrode positions. A mesh refinement may be needed, because a singular part would still exist in the case of a discontinuous conductivity model at an electrode location.

For a given topography, the Neumann conditions are not singular if the normals at the surface are computed between mesh nodes. The scalar product is indeed null close to the source. In the case where the normals are computed at the location of the potentials, that is, at node positions, the Neumann conditions are not singular is the curvature of the topography is small at the electrode positions (Appendix A).

Numerically, the secondary potential is indeed not singular (Fig. 1) for a locally homogeneous conductivity model and in the case where the topography is linearly interpolated between mesh points (blue curves) or when the curvature is small (red lines). A cubic interpolation is an intermediate solution for the description of the topography, but the curvature of the topography at the node locations and also in-between is higher than for the two other proposed solutions. This is why the secondary potential has larger values (green lines).

We conclude that a non-singular secondary potential can be obtained is the conductivity profile is locally homogeneous around the electrode positions and if the curvature of the topography is small at the position where the Neumann conditions are computed. This is not always the case, for example, if the conductivity profile has some discontinuities exactly at the source.

## IMPLEMENTATION

The approach proposed in Section 2 is valid for any numerical schemes (e.g. finite elements, finite differences, ...). We solve here eqs (11)–(13) with the GFDM (Liszka & Orkisz 1980). This method uses triangular and tetrahedral meshes (Geuzaine & Remacle 2009) to describe the conductivity model and the topography.

### Generalized finite difference method

GFDM generalizes finite difference formulae to define the derivatives of a function from a set of arbitrarily selected points. We consider here mesh connexions to define such a set (Fig. 2). For simplicity we present the method in 2-D but extension to 3-D is straightforward. For each point of a set localized around a central point (*x*_{c}, *z*_{c}), the Taylor series expansion of the potential *V* reads

*V*=

*V*(

*x*,

*z*),

*V*

_{c}=

*V*(

*x*

_{c},

*z*

_{c}),

*h*=

*x*−

*x*

_{c},

*k*=

*z*−

*z*

_{c}, $$\eta =\sqrt{h^{2}+k^{2}}$$. Eq. (14) should be valid for small η values and for all the

*m*points around (

*x*

_{c},

*z*

_{c}), including the central point. We thus obtain a set of

*m*linear equations

*x*

_{c},

*z*

_{c})

**A**contains all the geometrical parameters and is given by

The solution of this usually overdetermined set of linear equations is obtained by minimizing in the least-squares sens the norm ||**A****d _{v}** −

**v**||

^{2}, leading to

**A**

^{T}

**A**

**d**=

_{v}**A**

^{T}

**v**. We are less interested in the values of the partial derivatives and more in the weighting factors

*W*to apply to the points to obtain the finite difference approximation of the derivatives. The weights are defined by

**d**=

_{v}**W**

**v**, leading to

**W**= (

**A**

^{T}

**A**)

^{− 1}

**A**

^{T}, with

According to Liszka & Orkisz (1980) and Seibold (2006), we do not introduce the central point in the estimation of the derivatives for a better-conditioned matrix **A ^{T}A**. In practice, the stencil should contain the central point and at least five points in 2-D and nine points in 3-D to be well constrained.

### Conductivity averaging

Large conductivity contrasts, typically with ratios of 10 to 10^{4}, deteriorate the quality of the finite difference solution. To tackle this problem, we use a modified version of the Matched and Interface Boundary (MIB) method proposed by Yu *et al.* (2007), and combine it with the GFDM approach. In the case of unstructured grids, we obtain the classical formula as given below (Penz *et al.*2011). The GFDM discretization of the Laplacian of the potential (eq. 11) can be written as

*k*th connected point associated to the potential

*V*

_{s}. The key point is the arithmetic average of the conductivity values $$\bar{\sigma }_{k}$$ computed between two adjacent triangles or tetrahedrons. $$\bar{\sigma }_{k}$$ denotes the average conductivity of all cells containing the central point

*c*and a connected point

*k*.

### Free surface conditions

A specific treatment is applied for schemes with a central point located at the free surface. First for elements having an edge (or face) lying at the surface, we create fictitious symmetric triangles (or tetrahedrons) with respect to this edge (or face; Fig. 3) . Points 5 and 6 in Fig. 3 are considered to be fictitious points. The associated potentials are given by eq. (12). In the discrete domain, a null normal derivative reads *V*_{s5} = *V*_{s3} − (*V*_{p5} − *V*_{p3}) and *V*_{s6} = *V*_{s2} − (*V*_{p6} − *V*_{p2}). The primary potentials *V*_{p5} and *V*_{p6} are analytically known. There values only depends on the distance to the source. Finally, we replace *V*_{s5} and *V*_{s6} in eq. (20) to estimate ∇ · σ∇*V*_{s} at the central point.

### Solving the linear set of equations

Once the derivatives are locally estimated by minimizing *J*, eqs (11)–(12) are discretized, leading to a large sparse linear system. For *n* mesh nodes, the total size of the **A** matrix is *n* × *n*. However, for each line in the matrix, the number of non-zero elements is 1 + *m*, where *m* is the number of connected nodes to the central point. Typically, *m* varies between 5 and 10. The solution is obtained with a safe variant of the biconjugate-gradient method. A Jacobi pre-conditioner appears to help the convergence (Fig. 4).

## NUMERICAL EXAMPLES AND DISCUSSIONS

In this section, first we validate the proposed implementation by comparing the results for flat topography with analytical solutions. Then, we investigate the effect of topography on secondary potentials.

### A vertical contact model

The first model consists of a vertical contrast between a 100 Ωm medium and a 500 Ωm medium. Fig. 5 displays the apparent resistivity ρ_{a} of a Schlumberger configuration for a profile perpendicular to the contrast. The current electrode spacing is 24 m and the potential electrode spacing is 4 m. The 2.5-D and 3-D numerical solutions show satisfactory results: their curves nicely superimpose with the analytical solution. As expected, the 3-D code provides better results than the 2.5-D code. This can be explained by the limited number of wavenumbers (12) selected for the Fourier transform. The minimum wavenumber is set to 0, and the maximum wavenumber to 1/2/*dy* according to the Nyquist theorem, where *dy* is the lateral sampling that would be used in a 3-D setting (Xu *et al.*2000; Pidlisecky & Knight 2006). In our case, the wavenumbers are sampled with a third-order polynomial function, with more small wavenumbers. Additional wavenumbers would provide a better solution but at a higher computational cost, directly proportional to the number of selected wavenumbers.

### A buried sphere model

The apparent resistivity profile computed with the 3-D code in the case of a buried 0.5 Ωm conductive sphere in a homogeneous half space of 10 Ωm also shows a nice agreement with the analytical solution (Fig. 6). For this experiment, two fixed current electrodes are spaced by 600 m and two potential electrodes separated by 2 m are moved in-between and are never closer than 239 m from a current electrode.

### Topography effects

Figs 7– 10 display the values of the secondary potentials and of the total potentials for several source locations placed at the surface in an undulating topography case. The shape of the secondary potentials, computed in a homogeneous model (Figs 7 and 8), depends on the source location and on the topography. Most of the energy of *V*_{s} is concentrated around the source position due to the decay of the total potential as 1/*r*. Moreover, the sign of the secondary potential can be *a priori* estimated, but only on very simple cases (e.g., Fig. 7): a source close to a hill (respectively a valley) exhibits a positive (respectively negative) secondary potential. Note that with the *V*_{s} contribution, equipotentials of the total potential *V* are indeed perpendicular to the surface, as imposed by the Neumann conditions. The same conclusions hold for a 3-D modelling (Fig. 8). The rough aspect of the contour lines in the deeper part is due to a 3-D coarse grid discretization (Fig. 8, right-hand side).

For a void sphere buried in the subsurface, the structure of the secondary potential is much more complex, even for homogeneous conductivity models (Figs 9 and 10). For the same topography, but in the presence of a buried sphere, the secondary potential can be 3 to 4 times larger (Figs 9, left-hand side and 10, left-hand side). This happens when the source is close to the anomaly in the subsurface.

For a buried source, the same modelling strategy is applicable using *V*_{p} = *I*/(4πσ_{0}*r*) instead of *I*/(2πσ_{0}*r*) (Fig. 11). In this particular case, the secondary potential is positive everywhere.

A secondary potential for a source located close to a cavity is largely affected by the topography and by the presence of the cavity (Fig. 12, left-hand side). However, if the cavity is not at the vertical position of the electrode, the topography effect can be distinguished from the cavity effect (Fig. 12, right-hand side). The ratio *V*_{s}/*V*_{p} is similar to the case without cavity, except close to the cavity.

### Evaluation the topography effect

We propose two qualitative assessments of the topography effect. For simple topographies, we compute the sign of ∂_{n}*V*_{s} at the surface (Fig. 13). For an anticline structure and a source at the top of the anticline, the Neumann conditions are positive around the source position. This is the opposite for a syncline structure. For other source positions, we may have positive and negative values close to the source. These values can be checked by looking at the sign of the angle between the normal at the surface and the circular potential lines centred at the source position. This is a simple way to check the results obtained in Fig. 7, knowing that the key values are located around the source position.

A better approximation is provided with the fundamental solution method (FSM, Poullikkas *et al.*1998; Tsai *et al.*2006; Lin *et al.*2011). It consists of computing the secondary potential as a weighted sum of Green's functions for extra sources positioned outside of the zone of interest. The weights associated to the sources are defined such that the Dirichlet and Neumann conditions are satisfied at the border of the volume. The main advantage is that only the frontier needs to be discretized. The obtained solution is close to the numerical GFD solution (Fig. 14). The largest differences are observed at the vicinity of the source. The comparison however remains qualitative as it depends on the number and the position of extra sources (Lin *et al.*2011). Here, we selected three times more extra sources than collocation points used to described the boundary. We also impose smoothness constraints on the solution. A recommended solution would be to select the same number of sources and control points, but that solution is unstable with respect to the boundary conditions (Kirkup & Yazdani 2008).

### Discussions

The method presented in Section 2.2 was illustrated by a series of numerical tests. The modified singularity removal technique remains simple, even in the presence of topography. It does not require the numerical computation of the singular potential (Rücker *et al.*2006; Blome *et al.*2009) and is therefore advantageous in terms of computational cost. The idea of reformulating the Neumann conditions was already formulated by Lowry *et al.* (1989) despite they limited their study to a flat topography case. More recently Blome *et al.* (2009) decomposed the total potential into three terms $$V = V_p + V_s^1 + V_s^2$$. The classical primary potential *V*_{p} equals to −*I*/(2πσ_{0}*r*), whereas $$V_s^1$$ satisfies $$\Delta V_s^1=0$$ and $$\partial V_s^1 = -\partial V_p$$. Finally, $$V_s^2$$ is the solution of $$\nabla \cdot \sigma \nabla V_s^2 =-\nabla \cdot (\sigma -\sigma _0)\nabla (V_p+V_s^1)$$ and $$\partial _{n}V_s^2 = 0$$. It can be easily checked that $$V_s=V_s^1+V_s^2$$. The secondary field $$V_s^1$$ takes into account the topography effect and is computed numerically with the Boundary Element Method (BEM, Okabe 1981). The BEM approach certainly provides an accurate solution, but requires a fine description of the top interface as the Green's functions, containing a singularity, have to be integrated. This strongly limits the applicability of the BEM method. In practice, the Green's functions are decomposed differently for the near and far fields utilizing the fast multipole algorithm to obtain a satisfactory solution, but this largely complicates the approach (Hackbush & Nowack 1989; Blome *et al.*2009).

The alternative solution presented here is simpler as there is no need to define two secondary potentials and more importantly it does not require to integrate singular functions. The gradient of *V*_{p} is computed analytically, whereas in Blome *et al.* (2009), one needs also to estimate numerically $$\nabla V_s^1$$ for the derivation of $$V_s^2$$. In the modified approach, the boundary conditions can be imposed to be smooth if needed for a correct description of the topography. This has obvious advantages from a numerical point of view. Furthermore, the same strategy holds whatever the numerical scheme used for solving the equations. It could, for example, very well be used with a finite element approach, the key point being the implementation of given flux conditions at the free surface.

Note that the source position does not necessarily have to be defined on a mesh node: the proposed algorithm indeed requires to compute *V*_{p} and its gradient on the mesh, but these values are known analytically, whatever the source position. Then *V*_{s} is evaluated on the mesh. The flexibility with respect to the source position is an advantage, especially for finite difference schemes.

## CONCLUSIONS

We have proposed to extend the approach proposed by Lowry *et al.* (1989) in the presence of topography. We simply took as primary potential the analytic solution valid for a flat topography. For locally homogeneous conductivity models and when the topography presents a small curvature, the secondary potential does not contain any singularities. We modified accordingly the Neumann conditions for the secondary potential that become non-zero in the presence of topography and/or for heterogeneous models. The strategy is easy to be implemented, for surface and buried electrodes, even with possible cavities in the subsurface. Such benefits should be exploited in a future work in resistivity tomography, even if for a discontinuous conductivity profile exactly at the electrode positions, the secondary potential still remains singular.

The project was funded through the M.I.N.E.S Carnot research program.

## REFERENCES

### APPENDIX: SINGULARITY AT THE SOURCE AND FOR THE NEUMANN CONDITIONS

We examine the possible singularity in the secondary potential, arising either at the source or in the Neumann conditions in the presence of topography. The source term is given by eq. (11). For a smooth conductivity function around the source, it is proportional to *S* = ∇(σ − σ_{0}) · ∇(1/*r*), because 1/*r* is the Green's function of the Laplace equation, where *r* is the distance to the source location. A limited development of the conductivity in *r* is given by σ ≃ σ_{0} + α*r* + β*r*^{2} + γ*r*^{3}. A source term *S* without singularities imposes α and β to be equal to zero. It means in practice that the conductivity should be very smooth, possibly homogeneous, close to the source. If this is not the case, the mesh should be refined to remove the singularity in the source term to decrease the distance *r*. This is not always possible, in particular when the source is located at the discontinuity location of the conductivity or if the conductivity has a non-zero spatial gradient around the source.

For this particular case where (α and β ≠ 0), it is quite often that σ − σ_{0} ≠ 0 over subregions around the source point. Compared with the original total potential approach where the source is injected into a single geometrically singular point, the enlarged subregions over which σ − σ_{0} ≠ 0 on the secondary approach can significantly reduce the singular behaviour of the fields cause by the right-hand side source term in eq. (11).

The second cause for singularity is in the Neumann conditions around the source as indicated by eq. (12). The conditions are proportional to *N* = (**x** − **x _{s}**) ·

**n**/

*r*

^{3}. We examine the case when

*N*is not singular. We first consider that the normal and the potentials are computed at the same points. We apply a rotation centred around the source position such that the topography becomes locally horizontal. For a smooth topography and after rotation, the elevation

*z*=

*z*+ (

_{s}**x**−

**x**) ·

_{s}**n**is function of the position

*x*to the source

*x*

_{s}as

*z*≃

*z*

_{s}+ α(

*x*−

*x*

_{s}) + β(

*x*−

*x*

_{s})

^{2}+ γ(

*x*−

*x*

_{s})

^{3}. By definition α = 0 due to the rotation. The

*r*distance can be approximated by |

*x*−

*x*

_{s}|. The

*N*terms thus becomes |

*N*| ≃ β/

*r*+ γ. The Neumann conditions are not singular if β = 0, meaning that the curvature of the topography should be small at the electrode positions.

The second possibility is to compute the normals between two grid nodes. In that case, a linear interpolation between grid nodes at the surface ensures that the Neumann conditions are null around the source as the vector product is null.