2. Fundamental theory

According to the continuity of electrical current, the following integral equation is satisfied at any point in a conductive medium:

$$
\oint\_{\Gamma} \mathbf{J} \cdot d\Gamma = -I,\tag{1}
$$

where Γ is a full or half spherical surface that encloses an electrode that injects electric current J with a magnitude of I (Figure 1). According to Ohm's law:

$$
\mathbf{J} = \sigma \mathbf{E} = -\sigma \nabla U,\tag{2}
$$

and the property of the delta-function δ(x-xs), Eq. (1) gives

$$-\oint\_{\Gamma} \sigma \nabla U \cdot d\Gamma = -I \iiint\_{\Omega} \delta(\mathfrak{x} - \mathfrak{x}\_{i}) d\Omega \tag{3}$$

Here σ is the conductivity of medium, U is the electric potential, Ω represents the volume of the spherical surface Γ, and x<sup>s</sup> = (xs,ys,zs) is the location of the current electrode. Applying the divergence theorem to Eq. (3) leads to

Ω

ððð ½∇ � ðσ∇UÞ þ Iδðx�xsÞ� ¼ 0: (4)

Figure 1. Electric current density J around an electrode in (a) a full space and (b) a half-space.

Electrical Resistivity Tomography: A Subsurface-Imaging Technique DOI: http://dx.doi.org/10.5772/intechopen.81511

Note that Eq. (4) satisfies everywhere in a medium, so that the following governing equation of electric field is obtained:

$$\nabla.(\sigma \nabla U) = -I \mathfrak{K}(\mathfrak{x} - \mathfrak{x}\_{\mathfrak{t}}), \mathfrak{x} \in \Omega. \tag{5}$$

In general, the conductivity σ(x) = σ(x,y,z) changes with three coordinates so that it defines a 3D geological model, and the electric potential U=U(x,y,z) becomes a 3D function of a given current injection Iδðx � xsÞ and conductivity σ(x) (Figure 2a). If the conductivity is a constant with the y-coordinate, i.e. σ(x) = σ(x,z), it defines a 2D geological model (Figure 2b), but the electric potential U(x,y,z) is still a 3D function of the coordinates due to the point current injection Iδðx � xsÞ. Applying cosine Fourier transform Fc{�} to the y-coordinate of Eq. (5), the governing equation is changed into

$$\nabla \cdot \left( \sigma \nabla \tilde{\mathbf{U}} \right) + k\_\gamma^2 \tilde{\mathbf{U}} = -I \delta(\mathfrak{x} - \mathfrak{x}\_\circ)/2, \mathfrak{x} \in \Omega,\tag{6}$$

� � � � which is often named 2.5D governing equation. Here ky is the wavenumber, and Ue x; ky; z ¼ FcfUðx; y; zÞg becomes the spectrum of electric potential in the wavenumber domain. The coordinate vector x = (x,z) and gradient ∇ = (∂x, ∂z) are the 2D versions. If a geological model σ(x,y,z) or σ(x,z) and current injection <sup>I</sup>δð<sup>x</sup> � <sup>x</sup> <sup>Þ</sup> are given, solving the governing equation (5) or (6), one obtains electrical <sup>s</sup> � � potential U(x,y,z) for <sup>a</sup> 3D conductivity model or its spectrum <sup>U</sup><sup>e</sup> <sup>x</sup>; ky; <sup>z</sup> for <sup>a</sup> 2D conductivity model and then performs the inverse cosine Fourier transform U(x, y, z) n o <sup>=</sup> <sup>F</sup>�<sup>1</sup> <sup>U</sup><sup>e</sup> <sup>x</sup>; <sup>k</sup>y; <sup>z</sup> to obtain electric potential U(x,y,z). These computations are <sup>c</sup> called forward modelling. In theoretical computation or numerical forward modelling, Green's function G or Ge—the electric potential response to a unit current injection (I = 1 Am)—is often applied, so as to remove the magnitude of electric current I in Eqs. (5) and (6). One may focus on the computation of Green's function, and then the electric potential U or Ue is computed by the multiplication of Green's function G with the magnitude of practical electric current I, i.e. U = I�G.

The simplest geological model is a homogenous half-space. Applying Eq. (3) (equivalent to Eq. (5)) to a constant medium, the surface integral is calculated by

$$
\oint\_{\Gamma} \sigma \nabla U \cdot d\Gamma = \sigma \nabla U \left(2\pi r^2\right) = -I,\tag{7}
$$

which gives the electric potential at distance r from a current source:

Figure 2. Sketch of (a) a 2D and (b) a 3D geological model defined by conductivity σ(x,z) and σ(x,y,z), respectively.

Applied Geophysics with Case Studies on Environmental, Exploration and Engineering Geophysics

$$\mathbf{U} = \frac{I}{2\pi\sigma r} = \frac{\rho I}{2\pi r},\tag{8}$$

where ρ = 1/σ is resistivity. In practice, to inject electric current into the ground, a pair of current electrodes C1 and C2 must be employed; one is positive (+I) and another is negative (-I). Thus, the electric potential at a point P with a pair of current electrodes is calculated by.

$$\mathbf{U} = \frac{\rho I}{2\pi} \left( \frac{\mathbf{1}}{r\_{\rm C1P}} - \frac{\mathbf{1}}{r\_{\rm C2P}} \right). \tag{9}$$

Here, rC1<sup>P</sup> and rC2<sup>P</sup> are distances of the observed point P to two current electrodes C1 and C2. To measure the potential on the earth surface, a pair of potential electrodes P1 and P2 is also required (Figure 1). According to Eq. (9), one has the following expression for the electric potential difference between two potential electrodes:

$$
\Delta U = \frac{\rho I}{2\pi} \left[ \left( \frac{1}{r\_{\rm C1P1}} - \frac{1}{r\_{\rm C2P1}} \right) - \left( \frac{1}{r\_{\rm C1P2}} - \frac{1}{r\_{\rm C2P2}} \right) \right], \tag{10}
$$

from which apparent resistivity ρ<sup>a</sup> is defined as follows:

$$
\rho\_a = K \frac{\Delta \mathbf{U}}{I} = K \Delta \mathbf{G},
\tag{11}
$$

where ΔG = ΔU/I and K is called geometry factor of electrode array given by

$$K = 2\pi \left[ \left( \frac{1}{r\_{\rm C1P1}} - \frac{1}{r\_{\rm C2P1}} \right) - \left( \frac{1}{r\_{\rm C1P2}} - \frac{1}{r\_{\rm C2P2}} \right) \right]^{-1},\tag{12}$$

which depends on the positions of four electrodes. Different layouts of four electrodes have variable geometry factors and are often called electrode arrays. In the traditional electrode arrays, C2 and P2 may be set up very far away from C1 and P1,so that C2 and P2 are treated as remote electrodes in theory, and distances rC1P2, rC2<sup>P</sup>2, rC2<sup>P</sup>1, and rC2P<sup>2</sup> are supposed to be infinite (∞) in Eq. (12). In these cases, one can find that geometry factors are still applicable for these electrode arrays that involve one or two remote electrodes, which are named pole-pole, pole-dipole and dipole-pole arrays.

If subsurface resistivity is homogenous (ρ0), Eq. (11) shows that no matter which electrode array is used, apparent resistivity is constant (ρ<sup>a</sup> = ρ0). Otherwise, ρ<sup>a</sup> indicates resistivity variation of the underground. For a certain range of apparent resistivity ρa, Eq. (11) also reveals that the geometry factor is inversely proportional to the potential difference ΔG = ΔU/I. It implies that a big value of geometry factor K will cause a small reading of the electric potential difference ΔG in fields. Such a small reading is easily obscured or contaminated by background noise. For a good data quality, one should avoid very big values of geometry factors in data acquisition for ERT. Therefore, Eqs. (11) and (12) are the fundament formulae of the traditional DC electrical resistivity exploration.

## 3. Electrode arrays

In order to obtain apparent resistivity in fields, many electrode arrays were developed in the traditional DC resistivity exploration. In principle, ERT requires a Electrical Resistivity Tomography: A Subsurface-Imaging Technique DOI: http://dx.doi.org/10.5772/intechopen.81511

high data density and good coverage of the earth surface for high-resolution images of subsurface targets. Dahlin and Zhou [11] carried out synthetic experiments of ERT using 10 electric arrays and compared their imaging results for four geological models: a buried channel, a narrow dike, dipping blocks and waste ponds. They demonstrated that two three-electrode arrays (pole-dipole and dipole-pole) and three four-electrode arrays (dipole–dipole, Schlumberger and gradient arrays) produce satisfactory images of the subsurface targets. However, due to the use of remote electrodes, the three-electrode arrays are rarely applied for ERT in practice; thus, the four-electrode arrays become popular. Particularly, gradient array [12] is well suited for multichannel data acquisition and can significantly increase the speed of data acquisition in the field, and at the same time, it gives higher data density and lower sensitivity to noise than dipole-dipole array. Figure 3 shows three ˜ ° four-electrode arrays (upper row) and their pseudo-sections of data points <sup>x</sup><sup>ρ</sup><sup>a</sup> ; <sup>z</sup><sup>ρ</sup><sup>a</sup> (middle row) and their geometry factors (bottom row) for a layout of total 81 electrodes on the earth surface. Figure 3 shows that dipole–dipole array has the biggest range of geometry factor among these three-electrode arrays, and gradient array performs the smallest range, so that it is well suitable for data acquisition with a high data density and good data quality. Zhou and Dahlin [13] based on many sets of field data obtained with different electrode arrays and found that the measured potential errors depend on the reading of potential difference, ΔG = ΔU/I, and the measured datum in field may be expressed by

$$
\Delta \mathbf{G}^\* = \Delta \mathbf{G} \left[ \mathbf{1} + \mathbf{R} \left( \frac{a}{\Delta \mathbf{G}} \right)^\beta / \mathbf{100} \right], \tag{13}
$$

where <sup>Δ</sup>G\* denotes <sup>a</sup> noise-contaminated datum, <sup>R</sup> is <sup>a</sup> random number in (�1,1) and α and β are constants. Due to having a big range of geometry factor (see Figure 3), the dipole-dipole array is more sensitive to noise than Schlumberger and gradient arrays, so that one must pay much attention to the noise contamination when conducting dipole-dipole measurements [13]. According to the relationship between the geometry factor K and reading of potential difference ΔG, Eq. (13) may be applied to predications of noise-contaminated data for different electrode arrays [11].

Reviewing Figure 3, one may find that gradient array only uses the potential electrodes between C1 and C2 to measure the potential differences. It does not do the measurements at the potential electrodes outside of C1 and C2. It means that gradient array actually misses the outside data for every pair of current electrodes. To obtain a better data coverage and complete subsurface information, the outside

#### Figure 3.

Dipole-dipole, Schlumberger and gradient arrays (upper row) for ERT data acquisition and examples of their pseudo-sections of data points (x<sup>ρ</sup><sup>a</sup> , z<sup>ρ</sup><sup>a</sup> Þ (middle row) and geometry factors (bottom row).

potential electrodes of each pair of C1 and C2 should be employed in the gradient measurements. These additional data apparently complement to common gradient measurements and may improve ERT imaging. Accordingly, a new electrode array is naturally formed and called full-range gradient array shown in Figure 4a. The difference from the common gradient array is that the new electrode array uses not only the potential electrodes between C1 and C2 but also the outside potential electrodes of every pair of C1 and C2, if their geometry factors fall in a reasonable range for the data acquisition. Figure 4b and c shows the pseudo-section points and geometry factors for a layout of total 81 electrodes. Comparing with Figure 3, the full-range gradient array does improve the data coverage, and its geometry factors are controlled in a reasonable range.

If there are boreholes in a field, cross-hole ERT may be carried out to image the geological structure between the boreholes. Zhou and Greenhalgh [14] investigated all possible electrode arrays for cross-hole ERT data acquisition and found that the electrode arrays of pole-pole (A-M), pole-bipole (A-MN), bipole-pole (AB-M), and bipole-bipole (AM-BN) with their multi-spacing cross-hole profiling and scanning surveys are useful for cross-hole ERT. Here the capital letters A and B stand for two current electrodes, and M and N denote two potential electrodes. These cross-hole electrode arrays are shown in Figure 5 with their sensitivity functions in backgrounds [15]. They also found that the electrode arrays which have either both current electrodes or both potential electrodes in the same borehole, e.g. A-MN, AB-M and AB-MN, have a singularity problem in data acquisition (geometry factor goes to infinite so that apparent resistivity and pseudoreaction are not applicable), namely, zero readings of the potential or potential difference in cross-hole measurements, so that the potential data are easily obscured by background noise and their images are inferior to those from other

#### Figure 4.

Full-range gradient array (a) for ERT data acquisition and examples of its pseudo-section of data points ˜ ° <sup>x</sup><sup>ρ</sup><sup>a</sup> ; <sup>z</sup><sup>ρ</sup><sup>a</sup> (b) and geometry factors (c) with <sup>a</sup> layout of total <sup>81</sup> electrodes.

Electrical Resistivity Tomography: A Subsurface-Imaging Technique DOI: http://dx.doi.org/10.5772/intechopen.81511

#### Figure 5.

Electrode arrays for cross-hole ERT data acquisition. A and B stand for two current electrodes. M and N denote two potential electrodes. The background contours are the sensitivity functions of the electrode array.

cross-hole electrode arrays. The data having the singularity problem may be predicated by zero values of the inverse geometry factors, which should be avoided in cross-hole data acquisition. Therefore, A-M, AM-B, and AM-BN with multi-spaces are recommended for cross-hole ERT.

## 4. Numerical modelling

To compute theoretical electric potential U(x,y,z) in a complex 2D or 3D geological model, one has to solve the governing equation (6) or (5) using a numerical approach and a computer. For computational efficiency, an underground artificial boundary Γ<sup>0</sup> is required to truncate an infinite geological model, so that the boundary condition on Γ<sup>0</sup> must be introduced to the governing equations. Applying Eq. (7) to the artificial boundary Γ<sup>0</sup> and zero normal component (J�n = 0) of electric current density J on the earth surface, one may find that the electric potential U holds

$$
\sigma \nabla U \cdot \mathbf{n} + \nu U = \mathbf{0},
\tag{14}
$$

where n is the normal vector of the artificial boundary Γ<sup>0</sup> and ν is calculated by ν = (r�n)/r 2 . Eq. (14) is the 3D boundary condition on Γ0. Similarly, one may find the boundary condition for a 2D geological model [16]:

$$
\sigma \nabla \ddot{\mathbf{U}} \cdot \mathbf{n} + \nu \ddot{\mathbf{U}} = \mathbf{0}.\tag{15}
$$

Adding the artificial boundary conditions Eqs. (14) and (15) to Eq. (5) and (6), numerical modelling becomes to solve the following definite governing equation:

$$\text{3D}: \begin{cases} \nabla \cdot (\sigma \nabla G) = -\delta(\mathfrak{x} - \mathfrak{x}\_{\mathfrak{t}}), \mathfrak{x} \in \Omega, \\ \sigma \nabla G \cdot \mathbf{n} + \nu \mathbf{G} = \mathbf{0}, \ \mathfrak{x} \in \Gamma\_0 \end{cases} \tag{16}$$

or

$$\text{2D}: \begin{cases} \nabla \cdot \left( \sigma \nabla \tilde{\mathbf{G}} \right) + \mathbf{k}\_{\mathcal{V}}^{2} \tilde{\mathbf{G}} = -\delta(\mathbf{x} - \mathbf{x}\_{s})/2, \mathbf{x} \in \Omega, \\ \qquad \sigma \nabla \tilde{\mathbf{G}} \cdot \mathbf{n} + \nu \tilde{\mathbf{G}} = \mathbf{0}, \ \mathbf{x} \in \Gamma\_{0}. \end{cases} \tag{17}$$

The numerical approach to the definite governing equations is called numerical modelling. The most popular numerical approaches are finite-difference and finite-element methods. The former is simple and straightforward, due to directly applying a finite-difference formula to the derivatives of the gradient ∇, so that the definite governing equations at all points in the model domain Ω and on the artificial boundary Γ<sup>0</sup> are discretised and assembled into a linear equation system [16–20]. The latter converts the definite governing equation into an integral equation by weighting residual principle or variational principle [21] and then carries out numerical volume integration to produce a linear equation system [22–26]. Therefore, no matter which approach is applied, the numerical modelling becomes to solve a linear equation system. For an instant, the following paragraphs briefly show an advanced 2D version of the finite-element method, called Gaussian quadrature grid (GQG) approach, because of its advantages over finite-difference method and other schemes of finiteelement approaches, e.g. high accuracy and easy implementation of discretisation of a geological model having arbitrary free-surface topography. For a 3D numerical modelling, one may follow the 2D procedures.

Weighting residual principle is to calculate the following integral of the governing equation in Eq. (17):

$$\int\_{\mathfrak{U}} \mathbf{W} \Big[ \nabla \cdot \left( \sigma \nabla \tilde{\mathbf{G}} \right) + k\_{\mathcal{\gamma}}^2 \tilde{\mathbf{G}} + \delta(\mathfrak{x} - \mathfrak{x}\_{\mathfrak{t}}) / 2 \Big] d\mathfrak{Q},\tag{18}$$

where W is an arbitrary weighting function. Applying the divergence theorem to the above and then submitting the artificial boundary condition yields

$$\iint\_{\Omega} \left( \sigma \nabla W \cdot \nabla \tilde{\mathbf{G}} + \sigma k\_{\mathcal{Y}}^2 W \tilde{\mathbf{G}} \right) d\Omega + \int\_{\Gamma\_0} \nu W \tilde{\mathbf{G}} d\Gamma = W(\mathfrak{x}\_{\mathbf{s}})/2. \tag{19}$$

In order to calculate the integrals, the model domain Ω is divided into a set of the no-overlap subdomains {Ωe, e = 1,2,…,Ne} that matches the free-surface topography and subsurface interfaces. In each subdomain, Gaussian abscissae {(xα, zβ), α,β = 1,2,…,NG} are employed (see examples shown in Figure 6) to discretise the subdomain and then Lagrange interpolation:

$$\widetilde{\mathbf{G}} = \sum\_{i,j=1}^{N\_G} l\_i(\mathbf{x}), l\_j(\mathbf{z}) \widetilde{\mathbf{G}}\_{\neq}^{(\epsilon)}, \tag{20}$$

and Gaussian weights {wαβ} are applied to calculation of the submain integrals. Consequently, Eq. (19) becomes

$$\begin{split} \sum\_{\epsilon} \left\{ \sum\_{a,\emptyset,i,j} \omega\_{a\emptyset} \sigma\_{a\emptyset} \sigma\_{a\emptyset} \left[ l\_i'(\mathbf{x}\_a) \delta\_{j\emptyset} \partial\_{\mathbf{x}} \mathbf{W}\_{a\emptyset} + \delta\_{ia} l\_j'(\mathbf{z}\_\emptyset) \partial\_{\mathbf{z}} \mathbf{W}\_{a\emptyset} + \mathbf{k}\_\mathcal{p}^2 \mathbf{W}\_{a\emptyset} \delta\_{ia} \delta\_{j\emptyset} \right] \widetilde{\mathbf{G}}\_{a\emptyset}^{(\epsilon)} \\ + \sum\_{\mathbf{a}} \nu\_{\mathbf{a}} \mathbf{W}\_{\mathbf{a}} \widetilde{\mathbf{G}}\_{\mathbf{a}}^{(\epsilon)} \right\} &= \mathbf{W}(\mathbf{x}\_i) / 2. \end{split} \tag{21}$$

ð Þ<sup>e</sup> � � Here, <sup>G</sup><sup>e</sup> αβ are the discrete values of the electric potential spectra <sup>G</sup><sup>e</sup> <sup>x</sup>; <sup>k</sup>y; <sup>z</sup> in the subdomain Ωe. Choosing the Lagrange basis polynomials as the weighting functions, i.e. {Wpq ¼ lpð Þ x lqð Þz , ∀p, qg, Eq. (21) is changed into

Electrical Resistivity Tomography: A Subsurface-Imaging Technique DOI: http://dx.doi.org/10.5772/intechopen.81511

Figure 6.

2D and 3D Gaussian quadrature grids for numerical modelling: (a) a 2D model, 8 � 10 subdomains and 7 � 7 Gaussian abscissae in each subdomain and (b) a 3D model, 10 � 10 � 10 subdomains and 5 � 5 � 5 Gaussian abscissae in each subdomain.

$$\begin{split} \sum\_{\mathbf{c}} \left\{ \sum\_{\mathbf{u}, \boldsymbol{\theta}, \boldsymbol{i}, \boldsymbol{j}} \mu\_{\alpha\beta} \sigma\_{\alpha\beta} \left[ l\_{i}^{\prime}(\mathbf{x}\_{a}) l\_{p}^{\prime}(\mathbf{x}\_{a}) \boldsymbol{\delta}\_{\beta\beta} \boldsymbol{\delta}\_{q\beta} + l\_{j}^{\prime}(\mathbf{z}\_{\beta}) l\_{q}^{\prime}(\mathbf{z}\_{\beta}) \boldsymbol{\delta}\_{ia} \boldsymbol{\delta}\_{pa} + \mathbf{k}\_{\beta}^{2} \boldsymbol{\delta}\_{pa} \boldsymbol{\delta}\_{q\beta} \boldsymbol{\delta}\_{ia} \boldsymbol{\delta}\_{\beta\beta} \right] \widetilde{\mathbf{G}}\_{a\beta}^{(\epsilon)} \\ + \sum\_{\mathbf{r}} \nu\_{\mathbf{r}} l\_{p}(\mathbf{x}\_{\mathbf{r}}) l\_{q}(\mathbf{z}\_{\mathbf{r}}) \widetilde{\mathbf{G}}\_{\mathbf{r}}^{(\epsilon)} \right\} = l\_{p}(\mathbf{x}\_{\mathbf{s}}) l\_{q}(\mathbf{z}\_{\mathbf{s}}) / 2, \end{split} \tag{22}$$

which can be rewritten in the matrix form for all points (p,q) and (i,j) in Ω:

$$\mathbf{M}\ddot{\mathbf{G}} = \mathbf{b}\_{\text{o}} \tag{23}$$

� � ð Þe where M is a square symmetric matrix assembled by the coefficients of Ge αβ in Eq. (22). Ge and b<sup>s</sup> are two vectors: the former consists of the discrete values of Ge at all points of the Gaussian quadrature grid, and the latter is a zero vector except for the component of 1/2 at the electric current source located at xs. Eq. (22) shows that the matrices M and Ge depend on the wavenumber k<sup>y</sup> and conductivity model σαβ. The wavenumber can be predicated by the minimum and maximum spaces of electrodes [17]. The linear equation system may be efficiently solved for multiple current sources B = {b1,b2,…,bN} by the Cholesky decomposition [27]. After obtaining the spectra Ge , the theoretical electric potentials are obtained by inverse cosine Fourier transform <sup>U</sup>ðx; <sup>y</sup>; <sup>z</sup>Þ ¼ <sup>I</sup> � <sup>F</sup>�<sup>1</sup> <sup>G</sup><sup>e</sup> . <sup>c</sup>

From Eq. (22), one can find that the model conductivity σ(x,z) is discretised by σαβ at the Gaussian abscissae instead of finite elements. The subdomain integrals are calculated by the weights wαβ to the discrete integrands. The dense abscissae and variable σαβ may describe the details of a complex geological model and generate high accurate solutions, so that sizes of the subdomains are not necessarily small but shaped to match the free-surface and subsurface interfaces. Apparently, accuracy of the numerical modelling depends on the number of the Gaussian abscissae in the subdomains. The more abscissae are applied, the more accurate results are yielded but cost more computer time. For efficient and accurate modelling, the minimum electrode space of electrode array may be chosen for the size of the subdomain, to which five Gaussian abscissae is applied to produce satisfactory solutions of numerical modelling. As examples, Figure 7 shows the pseudo-sections of dipole– dipole, Schlumberger and gradient arrays for a layout of total 141 electrodes over an anticline model. These results show that similar pseudo-sections of apparent resistivity are obtained by using the three electrode arrays. From the above analysis, one

Figure 7.

Numerical modelling for (a) an anticline model and the apparent resistivity pseudo-sections of (b) dipoledipole, (c) gradient and (d) Schlumberger arrays. The discrete anticline model is also given in the pseudosections.

can see that the finite-element approach does not calculate the high-order derivatives for the governing equation, but any finite-difference scheme does. Therefore, finite-element approach is often called a 'weak'solution of the governing equation against a 'strong'solution obtained by finite-difference method.
