Geophysical Prospecting 32, 124-131, 1984.

GEOELECTRICAL M O D E L CALCULATIONS FOR

E. M U N D R Y " "

ABSTRACT
MUNDRY,E. 1984, Geoelectrical Model Calculations for Two-Dimensional Resistivity Dis-
tributions, Geophysical Prospecting 32, 1 2 4 1 31.
For the calculation of geoelectrical model curves for a two-dimensional resistivity dis-
tribution, the potential equation is transformed by means of a Fourier cosine transform into a
two-dimensional Helmholtz equation containing the separation parameter 1.
The numerical solution of this equation for different values of 1 for an irregular grid is
obtained using the method of finite differences combined with the method of overrelaxation.
The method by which derivatives are replaced by finite differences turned out to be very
important, especially for high resistivity contrasts. After testing several methods designed to
deal with any type of resistivity distribution, a method of discretization similar to that used
by Brewitt-Taylor and Weaver (1976) for magnetotelluric modeling for H polarization was
found the best.
Examples are given of model curves for Schlumberger soundings over a vertical fault
covered by overburden. The incorrect use of horizontal-layer models leads to erroneous
interpretations that are more complex than the real subsurface situations.

1. INTRODUCTION
Suitable modeling methods are essential for quantitative or semi-quantitative inter-
pretation of geoelectrical resistivity measurements conducted in areas of compli-
cated subsurface conditions, e.g. in exploration for ore deposits or in tectonically
disturbed areas.
For three-dimensional bodies embedded in a homogeneous half-space, the
boundary element method (Okabe 1981) is a fast method for calculating geo-
electrical model curves. If the resistivity distribution is more or less arbitrary, the
method of finite differences or the method of finite elements must be used. The

* Paper read at the 44th meeting of the European Association of Exploration Geophysicists,
Cannes, France, June 1982, revision received April 1983.
** Niedersachsisches Landesamt fur Bodenforschung, Stilleweg 2, Postfach 51 01 53, 3000
Hannover 51, German F.R.

124
 GEOELECTRICAL M O D E L CALCULATIONS 125

routine use of these methods for three-dimensional models does not seem to be
practical given the present state of computer techniques (Dey and Morrison 1979b,
Scriba 1981). However, they can be used for two-dimensional models, thus provid-
ing a considerable improvement in interpretation in the case of many practical
applications.
Real two-dimensional distributions of electrical potential, as treated by Mufti
(1976, 1978, 1980), occur when a line source parallel to the strike direction of a
two-dimensional geological structure is used, or in special geological situations, e.g.
in measurements made above a rotation-symmetrical salt dome. Point sources,
which are by far the most often used in geoelectrics, produce a three-dimensional
field even with a two-dimensional structure. However, the calculation of these fields
can be reduced to a two-dimensional problem.

2. M E T H O D
For a two-dimensional distribution of electrical conductivity ~ ( xz), , (where the
x-axis is horizontal and normal to the strike direction and the z-axis is directed
downwards) and a point source A(x, ,0,O) of current strength I at the earth's surface
z = 0, the equation for the potential V(x, y, z)
v . (OVV) = -Z6(X - xs)6(y)6(z),
(where 6 is the Dirac function) can be transformed into a Helmholtz equation

ax E) + 2
"(. (.$) - A20Y + 16(x - xs)6(x) = 0
using a Fourier cosine transformation (Dey and Morrison 1979a)

P(x,Z; A) =
r V(X,y, Z) COS Ay dy.

The solution P of this equation as a function of A under the corresponding bound-

V(X,y,
il"
Z) = - P(x, Z; A) COS Ay dy,

thus permitting the apparent resistivity to be calculated for any given electrode
array.
The Helmholtz equation was solved numerically for a number of values of 1by
the method of finite differences and by successive overrelaxation. The discretization
method as described by Mufti (1976, 1978, 1980) and Dey and Morrison (1979a)
yields sufficiently accurate model curves only for relatively small resistivity ratios
(up to about 2 : l), as concluded from comparison with analytical solutions and
from results obtained with the boundary element method.
The most practical procedure-even for larger resistivity ratios-proved to be
126 E. M U N D R Y

one that was used in a similar manner for two-dimensional magnetotelluric model
calculations by Brewitt-Taylor and Weaver (1976) for H polarization. The values for
the function V are approximated for each grid point (with indices i, j ) in the x- and
z-directions by a parabolic interpolation. In addition, the conductivities oi + l,z, + l,z
at the center of the area elements are defined at the beginning of the iteration
procedure, and weighted averages oi, j , and (dr~/dz)~, are calculated
j , (~?o/C?x)~,
according to the proportion of the area element with respect to the total area.
Details are given in a report by Mundry (1981).

3. E X A M P L EA: V E R T I C A L F A U L T
WITH OVERBURDEN

Figures 1 and 2 show practical examples of Schlumberger model curves for a 2-D
structure : soundings over a fault covered by overburden with configurations perpen-
dicular (fig. 1) or parallel (fig. 2) to the fault.
Using a fault with uniform cover as an example, fig. 3 shows how the conven-
tional interpretation of the data-assuming horizontal layering-yields an incorrect
model of the subsurface in the vicinity of the fault. If the center of the Schlumberger
array (assuming a configuration perpendicular to the fault) is above the better
conducting part of the subsurface, the model with the vertical boundary appears to

(a) 1 I I I I I l l I I 1

Fig. 1. Model curves for a Schlumberger configuration perpendicular to a vertical fault with
overburden. Three distances x from the sounding point to the fault were assumed (x/d = 0, 1,
2, where d is the thickness of cover in (a), (b), and (c), respectively).
GEOELECTRICAL M O D E L CALCULATIONS 127

\
- \

- x/d=l
‘L.
---____ 1-0.2- 0.2-

(b)
0.1 - I I I I l l I I I

-
xld = 2
.--
-- -- _ _ 1-0.2-0.2-

(C)
I I I I I I l l - I I I
128 E. M U N D R Y

Fig. 2. Model curves for a Schlumberger configuration parallel to a vertical fault with over-
burden. Three distances x from the sounding point to the fault were assumed ( x / d = 0, 1, 2,
where d is the thickness of cover in (a), (b), and (c), respectively).

(b)
 GEOELECTRICAL MODEL CALCULATIONS 129

(C)

be a three-layer case with a relatively low resistivity in the middle layer (H-type).
The thickness of this middle layer can be varied due to the S equivalence. The
relatively complicated interpretation is recognizable, although the data are deter-
mined with a relatively simple subsurface model.
For a configuration parallel to the fault, the model curves for sounding points 1
and 2 (above the low conductivity part) appear as K-type curves with some T
equivalence. The maximum deviations of the curves for the two models are of the
order of 1-3%.

ACKNOWLEDGMENTS
The work discussed in this paper was carried out within the scope of a research
project supported by the Federal Ministry for Research and Technology (Project
No. NTS 105 2).
Computer time was provided by the Regional Computer Center for Lower
Saxony in Hanover (Control Data Cyber 76).
The author would also like to express his appreciation to Dr P. Weidelt for
suggestions and discussions, Mr R. Schulz for providing a comparison with the
boundary element method, and Dr R.C. Newcomb for the translation.
130 E. M U N D R Y

1 2 3 4 5
0-
m
100

I0 -

20 -
-7-
30 -

100 100 100 100 100

v-
---lr---
370 340

perpendicular
mnm

1 2 3 4 5
0-
m
100 100 100 100 100

10 -

20 -

30 -

40 -

Fig. 3. Top: a model for a vertical fault with overburden. It is assumed that soundings with a
Schlumberger configuration are carried out at points 1-5. The corresponding sounding curves
can be taken from figs 1 and 2.
Center: interpretation of model curves for the model at the top by assuming a horizontal
layer case for each sounding point for a Schlumberger configuration perpendicular to the
fault.
Bottom: interpretation of model curves for the model at the top by assuming a horizon-
tal layer case for each sounding point for a Schlumberger configuration parallel to the fault.
 GEOELECTRICAL M O D E L CALCULATIONS 131

REFERENCES
BREWITT-TAYLOR, C.R. and WEAVER, J.T. 1976, On the finite difference solution of two-
dimensional induction problems, Geophysical Journal of the Royal Astronomical Society
47,375-396.
DEY, A. and MORRISON,H.F. 1979a, Resistivity modelling for arbitrarily shaped two-
dimensional structures, Geophysical Prospecting 27, 106-136.
DEY, A. and MORRISON, H.F. 1979b, Resistivity modelling for arbitrarily shaped three-
dimensional structures, Geophysics 44, 753-780.
MUFTI, I.R. 1976, Finite-difference resistivity modelling for arbitrarily shaped two-
dimensional structures, Geophysics 41, 62-78.
MUFTI,I.R. 1978, A practical approach to finite-difference resistivity modelling, Geophysics
43,930-942.
MUFTI,I.R. 1980, Finite-difference evaluation of apparent resistivity curves, Geophysical
Prospecting 28, 146-166.
MUNDRY,E. 1981, Geoelektrische Modellrechnungen fur zweidimensionale Strukturen,
Bericht BGR/NLfB 91 170, Hannover (unpublished report).
OKABE,M. 1981, Boundary element method for the arbitrary inhomogeneities problem in
electrical prospecting, Geophysical Prospecting 29, 39-59.
SCRIBA,H. 1981, Computation of the electric potential in three-dimensional structures, Geo-
physical Prospecting 29, 790-802.