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 xs = (xs, ys, zs) in a 3-D conductivity model σ(x). The electrical potential V(x) is governed by the Poisson equation and the following boundary conditions  

(1)
\begin{eqnarray} \nabla \cdot \sigma \nabla V &=& -I\delta (\mathbf {x}-\mathbf {x_{s}}), \end{eqnarray}
 
(2)
\begin{eqnarray} \partial _{n}V &=& 0 \rm {\quad at the ground surface}\ \Gamma , \end{eqnarray}
 
(3)
\begin{eqnarray} V &=& 0 \rm {\quad at the other boundaries}\ \partial \Omega , \end{eqnarray}
where δ is the Dirac delta function. The condition on the normal derivatives ∂nV 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 xs and a point in the subsurface x, they proposed the following mixed condition  
(4)
\begin{equation} \partial _{n}V+\frac{\cos \theta }{r} V=0 \rm {\quad at}\,\, \partial \Omega , \end{equation}
where θ is the angle between the radial vector and the outward normal to the boundary.

For a 3-D model invariant along the y-direction and a point source at (xs, 0, zs), 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  

(5)
\begin{equation} \tilde{V}(x,k_{y},z)=\int ^{\infty }_{0}V(\mathbf {x})\cos (k_{y}y)\,\mathrm{d}y, \end{equation}
where ky is the wavenumber related to y. The transformed potential is governed by the Helmholtz equation  
(6)
\begin{eqnarray} \nabla \cdot \sigma \nabla \tilde{V}(x,k_{y},z)-k_{y}^{2} \sigma \tilde{V}(x,k_{y},z)=-\frac{I}{2}\delta (x-x_{s})\delta (z-z_{s}). \nonumber\\ \end{eqnarray}
Eq. (6) has to be solved for different ky parameters. The 3-D potential is then obtained by performing the inverse Fourier transform  
(7)
\begin{equation} V(\mathbf {x}) = \frac{2}{\pi }\int ^{\infty }_{0}\tilde{V} (x,k_{y},z)\cos (k_{y}y)\,\mathrm{d}k_{y}. \end{equation}
The choice of the ky 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 Vp and a secondary potential Vs. The primary potential Vp contains the singularity due to the current source, whereas the secondary potential Vs contains the effects of a heterogeneous conductivity model. By definition, the total potential V is defined as V = Vp + Vs. Eqs (1) to (2) are rewritten as  

(8)
\begin{eqnarray} \nabla \cdot \sigma \nabla V_{s} &=& -\nabla \cdot \left[(\sigma - \sigma _{0})\nabla V_{p}\right], \end{eqnarray}
 
(9)
\begin{eqnarray} \partial _{n}V_s &=& 0 \rm {\quad at}\,\, \Gamma , \end{eqnarray}
 
(10)
\begin{eqnarray} V_s &=& -V_p \rm {\quad at}\,\, \partial \Omega , \end{eqnarray}
where Vp = I/2πσ0xxs‖ is the analytic solution valid for a homogeneous conductivity model and a source located at the surface. For mixed boundary conditions, eq. (10) becomes ∂nVs + cos θ/rVs = −∂nVp − cos θ/rVp. 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 ∂nVp = 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 ∂Vp = 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 Vp = I/2πσ0r remains identical. The only modification resides in the new definition of the Neumann conditions for the secondary potential  

(11)
\begin{eqnarray} \nabla \cdot \sigma \nabla V_{s} &=& -\nabla \cdot \left[(\sigma - \sigma _{0})\nabla V_{p}\right], \end{eqnarray}
 
(12)
\begin{eqnarray} \partial _{n}V_{s} &=& -\partial _{n}V_{p} = -\frac{I}{2\pi \sigma _0} \frac{(\mathbf {x-x_{s}}) \cdot \mathbf {n}}{r^{3}} \rm {\quad at}\,\, \Gamma , \end{eqnarray}
 
(13)
\begin{eqnarray} V_s &=& -V_p \rm {\quad at}\,\, \partial \Omega , \end{eqnarray}
where n is the outward normal at the surface. As above, for mixed boundary conditions, eq. (13) becomes ∂nVs + cos θ/rVs = −∂nVp − cos θ/rVp. By definition, ∂nVs=−∂nVp ensures that ∂nV = 0.

For flat topography, ∂nVp = 0 and thus eqs (9) and (12) are identical. For an arbitrary topography, Vp still contains the singular part but is analytically known. Vs 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).

Figure 1.

Shapes of the interpolated topography (green: after cubic interpolation, blue: after linear interpolation, red: after imposing a small curvature at the electrode positions) (top panel), and associated Neumann values ∂nVs for an electrode position at x = −8 m (bottom panel). The black points (top panel) indicate all fixed electrode positions.

Figure 1.

Shapes of the interpolated topography (green: after cubic interpolation, blue: after linear interpolation, red: after imposing a small curvature at the electrode positions) (top panel), and associated Neumann values ∂nVs for an electrode position at x = −8 m (bottom panel). The black points (top panel) indicate all fixed electrode positions.

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 (xc, zc), the Taylor series expansion of the potential V reads  

(14)
\begin{eqnarray} V - V_c &=& h\frac{\partial V_{c}}{\partial x}+ k\frac{\partial V_{c}}{\partial z}+ \frac{h^{2}}{2}\frac{\partial ^{2} V_{c}}{\partial x^{2}}+ \frac{k^{2}}{2}\frac{\partial ^{2} V_{c}}{\partial z^{2}} \nonumber \\ &&+\, hk\frac{\partial ^{2}V_{c}}{\partial x\partial z}+ O(\eta ^{3}), \end{eqnarray}
where V = V(x, z), Vc = V(xc, zc), h = x − xc, k = z − zc, $$\eta =\sqrt{h^{2}+k^{2}}$$. Eq. (14) should be valid for small η values and for all the m points around (xc, zc), including the central point. We thus obtain a set of m linear equations  
(15)
\begin{equation} \mathbf {A}\mathbf {d_v} = \mathbf {v}, \end{equation}
with  
(16)
\begin{equation} \mathbf {v}^{\rm T} =\lbrace V_{1},V_{2},\ldots ,V_{m} \rbrace , \end{equation}
and the unknown derivatives at the point (xc, zc)  
(17)
\begin{equation} \mathbf {d_v}^{\rm T} = \left\lbrace \frac{ \partial V_{c}}{\partial x}, \frac{\partial V_{c}}{\partial z}, \frac{\partial ^{2} V_{c}}{\partial x^{2}}, \frac{\partial ^{2} V_{c}}{\partial z^{2}}, \frac{\partial ^{2} V_{c}}{\partial x\partial z}\right\rbrace . \end{equation}
The matrix A contains all the geometrical parameters and is given by  
(18)
\begin{equation} \mathbf {A} = \left\lbrace \begin{array}{ccccc}h_{1} &\quad k_{1} &\quad\frac{h_{1}^{2}}{2}&\quad\frac{k_{1}^{2}}{2}&\quad h_{1}k_{1} \\ h_{2} &\quad \cdots &\quad \cdots &\quad \cdots &\quad\cdots \\ \vdots &\quad &\quad &\quad &\quad \\ \vdots &\quad &\quad &\quad &\quad \\ h_{m} &\quad &\quad &\quad &\quad \end{array} \right\rbrace . \end{equation}

Figure 2.

Selected points within a mesh for the estimation of the derivatives at the central point c.

Figure 2.

Selected points within a mesh for the estimation of the derivatives at the central point c.

The solution of this usually overdetermined set of linear equations is obtained by minimizing in the least-squares sens the norm ||Advv||2, leading to ATAdv = ATv. 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 dv = Wv, leading to W = (ATA)− 1AT, with  

(19)
\begin{equation} \mathbf {W} =\left[\begin{array}{lllll}w_{c}^{x}&\quad w_{1}^{x}&\quad w_{2}^{x}&\quad \cdots &\quad w_{m}^{x} \\ w_{c}^{z}&\quad w_{1}^{z}&\quad w_{2}^{z}&\quad \cdots &\quad w_{m}^{z}\\ w_{c}^{xx}&\quad w_{1}^{xx}&\quad w_{2}^{xx}&\quad \cdots &\quad w_{m}^{xx}\\ w_{c}^{zz}&\quad w_{1}^{zz}&\quad w_{2}^{zz}&\quad \cdots &\quad w_{m}^{zz}\\ w_{c}^{x}&\quad w_{1}^{xz}&\quad w_{2}^{xz}&\quad \cdots &\quad w_{m}^{xz}\\ \end{array}\right]. \end{equation}

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 ATA. 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 104, 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  

(20)
\begin{equation} \nabla \cdot \sigma \nabla V_{s} = -\left(\sum _{k=1}^{m}w_{k}\bar{\sigma }_{k} \right)V_{c}+ \sum _{k=2}^{m}w_{k}\bar{\sigma }_{k}V_{sk}, \end{equation}
where $$w_{k} = w^{xx}_{k}+w^{zz}_{k}$$ are the GFDM weights related to the kth connected point associated to the potential Vs. 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 Vs5 = Vs3 − (Vp5 − Vp3) and Vs6 = Vs2 − (Vp6 − Vp2). The primary potentials Vp5 and Vp6 are analytically known. There values only depends on the distance to the source. Finally, we replace Vs5 and Vs6 in eq. (20) to estimate ∇ · σ∇Vs at the central point.

Figure 3.

Selected points for the estimation of the derivatives when the central points belong to the free surface.

Figure 3.

Selected points for the estimation of the derivatives when the central points belong to the free surface.

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

Figure 4.

Convergence rates for different iterative solvers: Bi-Conjugate Gradient (BiCG), Conjudate-Gradient (CG) and a safe variant of the BiCG (BiCGSafe). The solid lines are without pre-conditioner and the dashed lines with a Jacobi pre-conditioner. The number of mesh nodes is 1386.

Figure 4.

Convergence rates for different iterative solvers: Bi-Conjugate Gradient (BiCG), Conjudate-Gradient (CG) and a safe variant of the BiCG (BiCGSafe). The solid lines are without pre-conditioner and the dashed lines with a Jacobi pre-conditioner. The number of mesh nodes is 1386.

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.

Figure 5.

Apparent resistivity curves for a Schlumberger array over a vertical contact located at x = 0 m (dashed line: 2.5-D computation, solid line: 3-D computation and dotted line: analytic solution).

Figure 5.

Apparent resistivity curves for a Schlumberger array over a vertical contact located at x = 0 m (dashed line: 2.5-D computation, solid line: 3-D computation and dotted line: analytic solution).

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.

Figure 6.

Apparent resistivity curves for a buried conductive sphere. The solid line corresponds to a 3-D computation and the dotted line denotes the analytic solution.

Figure 6.

Apparent resistivity curves for a buried conductive sphere. The solid line corresponds to a 3-D computation and the dotted line denotes the analytic solution.

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 Vs 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 Vs 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).

Figure 7.

Secondary potential in a linear scale (left-hand side) and total potential in a logarithm scale (right-hand side) for several source locations, with the 2.5-D code. The model is homogeneous at 100 Ωm.

Figure 7.

Secondary potential in a linear scale (left-hand side) and total potential in a logarithm scale (right-hand side) for several source locations, with the 2.5-D code. The model is homogeneous at 100 Ωm.

Figure 8.

Secondary potential in a linear scale (left-hand side) and total potential in a logarithm scale (right-hand side) for several source locations, with the 3-D code. The model is homogeneous at 100 Ωm.

Figure 8.

Secondary potential in a linear scale (left-hand side) and total potential in a logarithm scale (right-hand side) for several source locations, with the 3-D code. The model is homogeneous at 100 Ωm.

Figure 9.

Secondary potential in a linear scale (left-hand side) and total potential in a logarithm scale (right-hand side) for several source locations, with the 2.5-D code. The model is homogeneous at 100 Ωm.

Figure 9.

Secondary potential in a linear scale (left-hand side) and total potential in a logarithm scale (right-hand side) for several source locations, with the 2.5-D code. The model is homogeneous at 100 Ωm.

Figure 10.

Secondary potential in a linear scale (left-hand side) and total potential in a logarithm scale (right-hand side) for several source locations, with the 3-D code. The model is homogeneous at 100 Ωm.

Figure 10.

Secondary potential in a linear scale (left-hand side) and total potential in a logarithm scale (right-hand side) for several source locations, with the 3-D code. The model is homogeneous at 100 Ωm.

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 Vp = I/(4πσ0r) instead of I/(2πσ0r) (Fig. 11). In this particular case, the secondary potential is positive everywhere.

Figure 11.

Secondary potential in a linear scale (left-hand side) and total potential in a logarithm scale (right-hand side) for several buried source locations, with the 2.5-D code. The model is homogeneous at 100 Ωm.

Figure 11.

Secondary potential in a linear scale (left-hand side) and total potential in a logarithm scale (right-hand side) for several buried source locations, with the 2.5-D code. The model is homogeneous at 100 Ωm.

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 Vs/Vp is similar to the case without cavity, except close to the cavity.

Figure 12.

Ratios between the secondary potential Vs and the fixed primary potential Vp (top panel: from Fig. 7, lines 2 and 4, bottom panel: from Fig. 9, lines 2 and 4).

Figure 12.

Ratios between the secondary potential Vs and the fixed primary potential Vp (top panel: from Fig. 7, lines 2 and 4, bottom panel: from Fig. 9, lines 2 and 4).

Evaluation the topography effect

We propose two qualitative assessments of the topography effect. For simple topographies, we compute the sign of ∂nVs 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.

Figure 13.

Sign of the Neumann conditions for simple topography shapes and different source positions indicated by the green point. The red values denote a positive sign and the blue value a negative sign.

Figure 13.

Sign of the Neumann conditions for simple topography shapes and different source positions indicated by the green point. The red values denote a positive sign and the blue value a negative sign.

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

Figure 14.

Numerical solutions obtained with the GFD method (top panel), with the FSM approach (middle panel) and differences in absolute value (bottom panel). The secondary potential is expressed in mV. The source is positioned at xs = (8.0, 5.3) m.

Figure 14.

Numerical solutions obtained with the GFD method (top panel), with the FSM approach (middle panel) and differences in absolute value (bottom panel). The secondary potential is expressed in mV. The source is positioned at xs = (8.0, 5.3) m.

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 Vp equals to −I/(2πσ0r), 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 Vp 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 Vp and its gradient on the mesh, but these values are known analytically, whatever the source position. Then Vs 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

Blome
M.
Maurer
H.
Kersten
S.
Advances on 3D geoelectric forward solver techniques
Geophys. J. Int.
 , 
2009
, vol. 
176
 (pg. 
740
-
752
)
Coggon
J.H.
Electromagnetic and electrical modeling by the Finite Element Method
Geophysics
 , 
1971
, vol. 
36
 
1
(pg. 
132
-
155
)
Dey
A.
Morrisson
H.F.
Resistivity modeling for arbitrarily shaped three-dimensional structures
Geophysics
 , 
1979
, vol. 
44
 
4
(pg. 
753
-
780
)
Dieter
K.
Paterson
N.R.
Grant
F.S.
Ip and resistivity type curves for three-Dimensional Bodies
Geophysics
 , 
1969
, vol. 
34
 (pg. 
615
-
632
)
Geuzaine
C.
Remacle
J.-F.
Gmsh: a 3-D finite element mesh generator with built-in pre- and post-processing facilities
Int. J. Numer. Methods Eng.
 , 
2009
, vol. 
79
 
11
(pg. 
1309
-
1331
)
Günther
T.
Rücker
C.
Spitzer
K.
Three-dimensional modelling and inversion of dc resistivity data incorporating topography – II. inversion
Geophys. J. Int.
 , 
2006
, vol. 
166
 
2
(pg. 
506
-
517
)
Kirkup
S.
Yazdani
J.
A gentle introduction to the boundary element method in MATLAB/FREEMAT
Proceedings of the Conference on Mathematical Methods, Computational Techniques and Intelligent Systems
 , 
2008
Greece
Corfu
(pg. 
46
-
52
)
Lin
J.
Chen
W.
Wang
F.
A new investigation into regularization techniques for the method for fundamental solutions
Math. Comput. Simul.
 , 
2011
, vol. 
81
 (pg. 
1144
-
1152
)
Liszka
T.
Orkisz
J.
The finite difference method at arbitrary irregular grids and its application in applied mechanics
Comput. Struct.
 , 
1980
, vol. 
11
 
1–2
(pg. 
83
-
95
)
Lowry
T.
Allen
M.B.
Shive
P.N.
Singularity removal: a refinement of resistivity modeling techniques
Geophysics
 , 
1989
, vol. 
54
 
6
(pg. 
766
-
774
)
Mufti
I.R.
Finite-difference resistivity modeling for arbitrarily shaped two-dimensional structures
Geophysics
 , 
1976
, vol. 
41
 
1
(pg. 
62
-
78
)
Okabe
M.
Boundary element method formulations theory arbitrary inhomogeneities problem in electrical prospecting
Geophys. Prospect.
 , 
1981
, vol. 
29
 
1
(pg. 
39
-
59
)
Penz
S.
Chauris
H.
Donno
D.
Finite difference resistivity modeling on unstructured grids with large conductivity contrasts
Proceedings of the 17th European Meeting of Environmental and Engineering Geophysics, Near Surface 2011
 , 
2011
UK
Leicester
Pidlisecky
A.
Knight
R.
A MATLAB 2.5-D electrical resistivity modeling code
Comput. Geosci.
 , 
2008
, vol. 
43
 (pg. 
1645
-
1654
)
Poullikkas
A.
Karageorghis
A.
Georgious
G.
Methods of fundamental solutions for harmonic and biharmonic boundary value problems
Comput. Mech.
 , 
1998
, vol. 
21
 (pg. 
416
-
423
)
Queralt
P.
Pous
J.
Marcuello
A.
2-D resistivity modeling: an approach to arrays parallel to the strike direction
Geophysics
 , 
1991
, vol. 
56
 
7
(pg. 
941
-
950
)
Ren
Z.
Tang
J.
3D direct current resistivity modeling with unstructured mesh by adaptive Finite-Element Method
Geophysics
 , 
2010
, vol. 
75
 
1
(pg. 
H7
-
H17
)
Rücker
C.
Günther
T.
Spitzer
K.
Three-dimensional modelling and inversion of DC resistivity data incorporating topography – I. modelling
Geophys. J. Int.
 , 
2006
, vol. 
166
 (pg. 
495
-
505
)
Hackbush
W.
Nowack
Z.
On the fast matrix multiplication in the boundary element method by panel clustering
Numer. Math.
 , 
1989
, vol. 
54
 (pg. 
463
-
491
)
Seibold
B.
M-Matrices in meshless finite difference methods
PhD thesis
 , 
2006
Germany
University of Kaiserslautern
Spitzer
K.
A 3-D finite-difference algorithm for DC resistivity modelling using conjugate gradient methods
Geophys. J. Int.
 , 
1995
, vol. 
123
 
3
(pg. 
903
-
914
)
Tsai
C.C.
Lin
J.
Young
D.L.
Atluri
S.N.
Investigations on the accuracy and condition number for the method of fundamental solutions
Comput. Model. Eng. Sci.
 , 
2006
, vol. 
16
 
2
(pg. 
103
-
114
)
Yi
M.-J.
Kim
J.-H.
Song
Y.
Cho
S.-J.
Chung
S.-H.
Suh
J.-H.
Three-dimensional imaging of subsurface structures using resistivity data
Geophys. Prospect.
 , 
2001
, vol. 
49
 
4
(pg. 
483
-
497
)
Yu
S.
Zhou
Y.
Wei
G.
Matched interface and boundary (MIB) method for elliptic problems with sharp-edged interfaces
J. Comput. Phys.
 , 
2007
, vol. 
224
 
2
(pg. 
729
-
756
)
Zhou
B.
Greenhalgh
S.
Finite element three-dimensional direct current resistivity modelling: accuracy and efficiency considerations
Geophys. J. Int.
 , 
2001
, vol. 
145
 (pg. 
679
-
688
)
Xu
S.-Z.
Duan
B.-C.
Zhang
D.-H.
Selection of the wavenumbers k using an optimization method for the inverse Fourier transform in 2.5D electrical modelling
Geophys. Prospect.
 , 
2000
, vol. 
48
 
8
(pg. 
789
-
796
)

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 + βr2 + γr3. 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 = (xxs) · n/r3. 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 = zs + (xxs) · n is function of the position x to the source xs as z ≃ zs + α(x − xs) + β(x − xs)2 + γ(x − xs)3. By definition α = 0 due to the rotation. The r distance can be approximated by |x − xs|. 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.