Automatic Interpretation of 

Schlumberger Sounding Curves, 

Using Modified Dar Zarrouk 

Functions

GJEOLOGICA ti SURVEY B It t L E T I N 131J-E

Automatic Interpretation of 

Schlumberger Sounding Curves, 

Using Modified Dar Zarrouk 

Functions

By ADEL A. R. ZOHDY

NEW TECHNIQUES IN D I R E CT - CUR R ENT 

RESISTIVITY EXPLORATION

GEOLOGICAL SURVEY BULLETIN 1313-E

Formulas defining two types of modified Dar ^arrouk curves 

are used to invert Schlumberger sounding curves, using an 

iterative procedure. The number of layers, which is equal to 

the number of points on the inverted curve, is reduced by 

automatically smoothing the corresponding D£ curve

UNITED STATES GOVERNMENT PRINTING OFFICE, WASHINGTON : 1975
UNITED STATES DEPARTMENT OF THE INTERIOR 

JAMES G. WATT, Secretary

GEOLOGICAL SURVEY 

Dallas L, Peck, Director

First printing 1975 

Second printing 1983

For sale by the Distribution Branch, U.S. Geological Survey, 

604 South Pickett Street, Alexandria, VA 22304
CONTENTS

	Page 

Abstract __________________________________________ El

Introduction _______________________________________ 2

Earlier developments __________________________________ 3

Outline of method ____________________________________ 6

Construction and inversion of DZ curves ______________________ 7

Modified Dar Zarrouk curves of the L type _____________________ 13

Inversion of MDZ-L curves ______________________________ 14

Modified Dar Zarrouk curves of the T type _____________________ 18

Inversion of MDZ-T curves ______________________________ 20

Inversion of VES curves ________________________________ 21

Inversion of incomplete and (or) distorted VES curves ______________ 26

Automatic smoothing and inversion of DZ curves _________________ 29

Field example ______________________________________ 33

Advantages and limitation of the method ______________________ 36

References cited _____________________________________ 37

ILLUSTRATIONS

Page

FIGURE 1. Two-layer DZ curves _________________________ E8

2. Schematic diagram for construction and inversion of DZ curves 9

3. Schlumberger VES curves and their corresponding DZ curves _ 10

4. Graph showing inversion of an A-type three-layer VES curve _ 11

5. Three-layer VES, DZ, and Hummel curves ____________ 14

6. MDZ-L curves _____________________________ 15

7. MDZ-T curves _____________________________ 19

8. Curves for five iterations in the inversion of an H-type VES curve 24

9. DZ curves for three- and nine-layer models ____________ 25

10-13. Graphs showing:

10. Inversion of an incomplete and distorted VES curve _ _ 28

11. Variation of fitting tolerance as a function of DZ slope 33

12. Automatic interpretation of a field VES curve _____ 34

13. Manual smoothing and inversion of detailed DZ curve

for layering shown in figure 12 ____________ 35

III

NEW TECHNIQUES IN DIRECT-CURRENT RESISTIVITY 

EXPLORATION

AUTOMATIC INTERPRETATION OF

SCHLUMBERGER SOUNDING CURVES,

USING MODIFIED DAR ZARROUK

FUNCTIONS

By ADEL A. R. ZOHDY

ABSTRACT

For horizontally stratified media, Schlumberger VES (vertical electrical sounding) 

curves with slopes of greater than -1 not only resemble but commonly almost coincide 

with their corresponding DZ (Dar Zarrouk) curves. By considering an n-point 

Schlumberger VES curve to be an n-layer DZ curve, and by solving for the layering from 

the DZ curve, we obtain a first approximation to the actual layering in the form of an nlayer section. However, the minimum slope for DZ curves is 1, whereas on some VES 

curves negative slopes may be as low as 3; therefore, these VES curves cannot be considered to be a first approximation to their corresponding DZ curves. Formulas are obtained for calculating two types of MDZ (modified DZ) curves whose positive and 

negative slopes are not limited to +1 and 1, respectively. Therefore, for any VES curve, 

there exists a corresponding MDZ curve that lies close to it. The automatic interpretation 

is made by an iterative procedure in which, for the first approximation, the observed VES 

curve is assumed to be the sought MDZ curve. This MDZ curve is solved for layer 

thicknesses and resistivities; then, by means of a convolution technique, a VES curve is 

calculated for the obtained layering, and the calculated VES curve is compared with the 

observed one. A second approximation to the sought MDZ curve is obtained by utilizing 

the differences between the observed and calculated VES curves. The iteration is continued until a match, within a prescribed fitting tolerance, is obtained between observed 

and calculated VES curves. The number of layers in the resulting model (detailed 

solution) is always equal to the number of points used to define the observed VES curve. 

Equivalent solutions composed of a fewer number of layers are determined by 

automatically smoothing the DZ curve of the detailed solution and inverting it. Special 

equivalent solutions that are subject to certain geologic or geoelectric constraints can be 

found by manual adjustments of the detailed n-layer DZ curve. Excellent automatic 

matches were obtained when the method was tested with several theoretical VES curves 

and with several hundreds of field VES curves of different forms. The average processing 

time per VES curve (extending over more than three logarithmic cycles) on the IBM 

360/65 computer is about 8 seconds. For distorted and incomplete field curves, the 

average processing time is approximately doubled.

El
E2 TECHNIQUES IN DIRECT-CURRENT RESISTIVITY EXPLORATION

INTRODUCTION

The interpretation of VES (vertical electrical sounding) curves is 

usually made by curve-matching procedures in which albums of 

theoretical curves (Compagnie Generale de Geophysique, 1963; Flathe, 

1963; Orellana and Mooney, 1966; Rijkwaterstaat, 1969) are used in 

conjunction with the auxiliary-point method of partial curve matching 

(Kalenov, 1957; Orellana and Mooney, 1966; Zohdy, 1965). Although 

this method is very educational and significant, it requires a great deal 

of practice, and at times it can be frustrating to the inexperienced interpreter. In fact, for geoelectric sections containing layers with small 

effective-relative resistances or small effective-relative conductances 

(Flathe, 1963; Zohdy and Jackson, 1973; Zohdy, 1974a), even the experienced interpreter may find it difficult to obtain satisfactory curve 

matches when he uses the auxiliary-point method. The use of DZ (Dar 

Zarrouk) curves (Maillet, 1947; Orellana, 1963) is valuable in improving the theoretical fit to an observed VES curve (Zohdy, 1974a), especially when the interpretation must be modified so that the interpreted geoelectric section would contain a different number of layers 

than the one initially assumed.

The direct interpretation, or inversion, of VES curves into layer 

thicknesses and resistivities has attracted the attention of geophysicists 

for the past 40 years. In late November 1972 the author presented a 

paper on automatic interpretation of Schlumberger sounding curves, 

using modified Dar Zarrouk functions, at the 42d annual international 

meeting of the Society of Exploration Geophysicists in Los Angeles, 

Calif. The present report is based on that talk and describes a fast 

method for the automatic interpretation of Schlumberger-type VES 

curves which are obtained over horizontally stratified, laterally 

homogeneous media. The method is intended for the processing of large 

numbers of VES curves on a digital computer in very short periods of 

time so that the interpreter may have at his disposal a basic solution, 

and simplified solutions, for each observed VES curve. The basic solution also may be modified manually, through the use of DZ curves, to 

obtain electrically equivalent solutions that will fit certain constraints 

on the layering in the survey area. The layering provided by the computer, moreover, may be considered as an interpretation made by an independent interpreter whose computations do fit the observed data 

and, therefore, often may point out a certain layering distribution 

which the interpreter may have overlooked.

The method was designed so that positive, realistic, layer resistivities 

and thicknesses are always obtained, even if the VES curve is distorted 

by lateral heterogeneities. For these distorted VES curves, although a 

practically perfect match between observed and calculated curves may 

be impossible (because the solution is based on a horizontally stratified,
INTERPRETATION OF SCHLUMBERGER SOUNDING CURVES E3

laterally homogeneous earth model), the calculated curve will fit the 

observed points as closely as possible, except for those points near the 

distorted segment.

EARLIER DEVELOPMENTS

Several methods have been proposed for the direct interpretation, or 

inversion of VES curves obtained over horizontally stratified, laterally 

homogenous media. Stefanesco and others (1930) derived the integral 

expressions

/ 

V(r)= ^ \± +2 \Q(\,k,d)J0(\r)dx] (1) Mr Jo J

=pi | l+2r2

for expressing the potential V(r) and the Schlumberger apparent resistivity jos(r) as a function of the Schlumberger electrode spac-, 

ing, r=AB/2; the first-layer resistivity, pi; the Bessel functions, Jo and 

J\\ and the Stefanesco kernel function, Q(X,k,d), which is expressed in 

terms of the layer depths, d\ the resistivity reflection coefficients, k, and 

the integration variable, X. If the TKF (total kernel function), T(X,d,k), 

also known as the Schlichter kernel (Vozoff, 1958), is defined by

; (3) 

then,

e(X,d,fe) = [(T(X,d,fe)/pi)-l]/2. (4)

Substituting equation 4 in equation 1 we get

r°° r °°

V(r) =&? \ i + I T(X' d' k)J0 (\r)d\- I J0 (Xr)rfx] . (5) 

27r I 1 pl 1 J

But^ according to the Weber-Lipschitz identity (Watson, 1962), r°°

£ = I «/o(Xr)dX, (6)

 'o 

and therefore (Sunde, 1949) equation 5 reduces to

/»°° 

Vfr)=jL I T(\.d.k)Jn (\r)d\. (7)
E4 TECHNIQUES IN DIRECT-CURRENT RESISTIVITY EXPLORATION

Differentiating equation 7 with respect to r and substituting in the 

equation

a(r)_-&r'SJy. (8)

we get

T(\,d,k) MXr) - XdX. 1'0

With an appropriate change in the variables r and X, the integral in 

equation 9 can be transformed into a convolution integral (Kunetz, 

1966).

Inasmuch as the kernel function, Tor 0, contains all the information 

that is necessary to define the geoelectric section (layer thicknesses and 

resistivities), Schlichter (1933) suggested the following procedure for 

the direct interpretation of resistivity sounding data:

1. Determine the kernel-function curve from the apparent resistivity 

curve,

2. Solve for the layering from the kernel function. 

The value of T(\,d,k) can be expressed in terms of the Schlumberger 

apparent resistivity by applying the inverse Hankel transformation of 

the Fourier-Bessel integral (Watson, 1962) to equation 9 so that

T(\,d,k) = p s(r))dr. (10)

f

 '

The numerical evaluation of T(\,d,k) from equation 10 has been 

studied by several investigators. Strakhov (1966a, 1966b, 1968) derived 

expressions for the calculation of T(\,d,k), and Strakhov and Karelina 

(1969) published coefficients for the transformation of a Schlumberger 

VES curve into its corresponding TKF curve, using convolution. In 

order to make this transformation, using the Strakhov-Karelina coefficients, the VES curve must be digitized at the logarithmic electrodespacing interval of (AB/2)i+l/(AB/2)i=e(ln 10/5) =*1.58, which is 

equivalent to five logarithmically equally spaced points per logarithmic 

cycle. The digitized apparent resistivities are convolved with the 

Strakhov-Karelina coefficients to obtain the total-kernel-function 

curve.

Ghosh (197 la) used linear-filter theory to obtain other sets of coefficients for transforming Schlumberger and Wenner sounding curves 

into their corresponding TKF curves. Ghosh's coefficients are 

applicable for the logarithmic electrode-spacing interval of 

of (AB/2)i+l/(AB/2)i=e(ln 10/3)-2.15, but, for a greater definition 

of the calculated TKF curve, his coefficients may be convolved with
INTERPRETATION OF SCHLUMBERGER SOUNDING CURVES E5

two sets of digitized points on the VES curve, which are separated by an 

interval of 1.46, to obtain points on the TKF curve which are separated 

by a factor 1.46 instead of 2.15.

The use of convolution for the calculation of the TKF curve from an 

observed VES curve requires the extrapolation of the VES curve to the 

left of the first and to the right of the last measured points. The values 

of the ordinates of the extrapolated VES points not only will affect the 

values of the calculated TKF curve at the TKF points corresponding to 

the first and last VES points but also will affect the values of several 

successive TKF points to the right and to the left of the first and last 

calculated points, respectively. Furthermore, a distortion over a small 

segment of the VES curve will cause a distortion over a correspondingly 

larger segment on the transformed TKF curve.

Koefoed (1965, 1966, 1968) used the raised kernel function H(\) instead of the Stefanesco kernel function 9(A) and subsequently 

(Koefoed, 1970) introduced graphs to accelerate the evaluation of the 

layer thicknesses and resistivities from the raised kernel function H(X). 

Other students of the direct interpretation of resistivity data in the 

kernel domain include Crous (1971), Meinardus (1967, 1970), Onodera 

(1960), Pekeris (1940), and Vozoff (1958).

One of the most attractive features of transforming a VES curve into 

its corresponding TKF curve (and finding the appropriate horizontally 

layered model for which a calculated TKF fits the transformed curve) is 

the ease with which a TKF curve can be calculated using Sunde's, 

Flathe's, or Vanyan's recursion formulas (Sunde, 1949; Flathe, 1955; 

Vanyan and others, 1962; Kunetz, 1966; Meinardus, 1967). Russian 

geophysicists have compiled albums of theoretical TKF curves and 

prepared equivalence nomograms similar to the Pylaev nomograms 

(Bhattacharya and Patra, 1968) for the interpretation of transformed 

TKF curves by the method of curve matching (Shkabarnia and 

Gritsenko, 1971).

Kunetz and Recroi (1970) pointed out some of the limitations of the 

interpretation of VES curves in the kernel domain, such as the risk of 

rapid amplification of errors in the initial data and the possibility, with 

certain techniques, of obtaining solutions that are physically absurd 

(for example, negative resistivities and (or) negative thicknesses). 

Furthermore, because the maximum negative slope for a TKF curve is 

 1, similar to DZ curves, the resolving power of TKF curves for Q-type 

(pi>p2>ps) and certain HK-type (pi>p2<ps>p4) sections is lower than 

the resolving power of the corresponding VES curves; therefore, the 

range of the principle of equivalence, in the kernal domain, is larger 

than it is for VES curves (Strakhov and Karelina, 1969). In other words, 

if two TKF curves seem to coincide (where one curve was obtained by 

the transformation of an observed VES curve and the other was
E6 TECHNIQUES IN DIRECT-CURRENT RESISTIVITY EXPLORATION

calculated for an assumed model), it does not mean, necessarily, that 

their corresponding VES curves will also almost coincide. Therefore, for 

interpretations in the kernel domain, the calculation of the VES curve 

for comparison with the observed VES curve is advisable, especially for 

VES curves with steeply descending branches. In this regard, although 

Kunetz and Recroi based their analysis on the determination and solution of a "nucleus" kernel function, they were the first to show the applicability of their technique by comparing reference or observed VES 

curves with calculated VES curves, and not by comparing kernel function curves.

Kunetz (1966) indicated that it is possible to calculate VES curves by 

convolving the TKF curve with a set of coefficients which he did not 

publish. Ghosh (1971b) published two sets of inverse filter coefficients 

for the calculation of Schlumberger and Wenner VES curves by the 

method of convolution. The convolution technique is the fastest known 

method for the calculation of a VES curve. For example, a VES curve of 

the Schlumberger type, which is defined by 18 points that are distributed over three logarithmic cycles and which corresponds to a section composed of about 20 layers, can be calculated using the Ghosh 

coefficients in less than about 0.5 second on the IBM 360/65 computer 

(Zohdy, 1974b). With the advent of this technique, one of the most attractive features of interpreting VES curves in the kernel domain, 

namely the speed with which a TKF curve can be calculated, was 

rendered less significant.

The use of numerical integration (Meinardus, 1967,1970) yields more 

accurate but slower (by several seconds) computations of VES curves 

than the use of convolution. The numerical integration method, 

however, is at least 10-20 times faster than the methods that are based 

on the summation of a slowly convergent series (Mooney and others, 

1966). Furthermore, both the convolution method and the numerical integration method are not dependent on the condition that the layers in 

the section must be of commensurate thicknesses (thicknesses that are 

whole multiples of a given unit of length).

OUTLINE OF METHOD

The automatic interpretation method presented here is based on inverting the VES curve without first transforming it into its corresponding TKF curve. Furthermore, unlike most of the direct interpretation methods, the interpreter does not have to make an initial 

assumption on the number of layers or about their resistivities and 

thicknesses. Instead, the VES curve is digitized (at the rate of six points 

per logarithmic cycle), and the number of layers is automatically fixed 

as equal to the number of points on the digitized curve.

For the first approximation, the points on the digitized VES curve are
INTERPRETATION OF SCHLUMBERGER SOUNDING CURVES E7

considered to be points on a modified Dar Zarrouk (MDZ) curve, which 

is solved for layer thicknesses and resistivities. The TKF curve for this 

layering is calculated by Sunde's recursion formula, and the VES curve 

is calculated by convolution, using Ghosh's coefficients. The calculated 

and observed VES curves are compared, and, through an iterative formula, a new MDZ curve is calculated, solved for layering, and a second 

VES curve is calculated. The iteration continues until a fit, within a 

prescribed tolerance, is obtained between observed and calculated VES 

curves. In this procedure the calculation of the VES curve by convolution requires the calculation of the TKF curve at (I/A) values that extend beyond the range of the corresponding (AB/2) values. However, 

because these TKF values are calculated for a given model by a recursion formula and are not calculated by the transformation of a VES 

curve, they are all accurately determined over the entire range of (I/A) 

values.

CONSTRUCTION AND INVERSION OF DZ CURVES

DZ curves for n-layer sections are calculated from the parametric 

equations

Pm .= 

and

where pm . and Lm . are the DZ resistivity and DZ depth; p; and HI are the

resistivity and thickness of the ith layer; and TI and Si are the 

transverse resistance and longitudinal conductance of the ith layer, 

respectively.

Figure 1 shows a set of two-layer DZ curves calculated for various 

values of M=p2/pi. It is important to note that the maximum positive 

and negative slopes for p2=°° and p2=0 are +1 and 1, respectively. 

Graphical methods were devised by Orellana (1963) and by Zohdy 

(1974a), for the construction of multilayer DZ curves.

The inversion of DZ curves into layer thicknesses and resistivities is 

made using the formulas

'"V ' ~"V-i^-i (13)

L/m- Li m . ,

Pmj
E8 TECHNIQUES IN DIRECT-CURRENT RESISTIVITY EXPLORATION

100h0.01 I 

0.004 I I I I I I 

0.5 1 500

FIGURE 1. Two-layer DZ curves for various values of the resistivity ratio p2 /pi. Dashed 

curves designate the thickness ratio hi/hi, pm, DZ resistivity; Lm, DZ depth; pi and hi, 

first-layer resistivity and thickness.
INTERPRETATION OF SCHLUMBERGER SOUNDING CURVES E9

Pm;-m/_ 1 'S-l (14)

The above relationships are schematically summarized in figure 2; 

practical applications of the inversion of DZ curves (pertaining to the 

interpretation of VES curves) were given by Zohdy (1974a).

Figure 3 shows that Schlumberger VES curves on which there are no 

negative slopes of less than or equal to - 1 commonly almost coincide 

with their corresponding DZ curves. If these DZ curves are solved for 

layering, using equations 13 and 14, either by using only the terminal 

DZ points (which mark the end of one layer and the beginning of 

another) or by digitizing the DZ curves so that the terminal points are 

included, we obtain the exact geoelectric sections for which both the DZ 

and VES curves were calculated. Therefore, because of the near coincidence of these VES curves with their corresponding DZ curves, a first 

approximation of the actual layering can be obtained by digitizing the 

VES curves, and substituting, in the inversion equations 13 and 14, the

10,000

FIGURE 2. Schematic diagram of the construction and inversion of DZ curves, hj, layer 

thickness; PJ, layer resistivity; T, transverse resistance; S, longitudinal conductance; 

pm, DZ resistivity; Lm, DZ depth.
E10 TECHNIQUES IN DIRECT-CURRENT RESISTIVITY EXPLORATION

1000

.100

o 5

CL ul

UJ CC 10

CD§s

-J Q 

I Z

EXPLANATION

3 10 100 1000 2000 

SCHLUMBERGER ELECTRODE SPACING (AB/2) AND DZ DEPTH (Lm ), IN METRES

FIGURE 3. Graphs of theoretical Schlumberger VES curves and their corresponding DZ

curves.

values of the electrode spacings, AB/2, for the DZ depths, Lm, and the 

Schlumberger apparent resistivities, ps, for the DZ resistivities, pm, 

respectively.

Figure 4 shows five steps for the successive approximation of a layering for a theoretical VES curve. Step 1 shows the layering, the VES 

curve, and the DZ curve for an A-type section (pi<p2<pa). By digitizing 

the VES curve, as shown in step 2, assuming that the VES curve is a DZ 

curve, and solving for the layering, we obtain the rc-layer section shown 

by the dashed line. This layering serves as a first approximation to the 

true layering. By calculating the TKF curve for the n-layer section from 

Sunde's recursion formula and convolving it with Ghosh's coefficients, 

the VES curve for this layering can be calculated and plotted as shown 

by the dashed curve of VES' in step 3. By knowing the effect of changes
INTERPRETATION OF SCHLUMBERGER SOUNDING CURVES Ell

1000p i i i i i n

10 100 1000 

AB/2,Lm, DEPTH, IN METRES

10,000

FIGURE 4. Graphs showing the inversion of a three-layer VES curve of the A-type by 

iteration through successive approximations of its corresponding DZ curve. Circles on 

VES curve designate first approximation to DZ curve. Arrow heads designate second 

approximation to DZ curve.
E12 TECHNIQUES IN DIRECT-CURRENT RESISTIVITY EXPLORATION

in the form of DZ curves on the form of their corresponding VES curves, 

as they were recently described by Zohdy (1974a) we could have 

predicted the approximate shape of the VES' curve in relation to the 

shape of the reference VES curve by noting the deviations of the ordinates of the reference VES curve from the ordinates of its corresponding DZ curve.

In order to obtain a second approximation to the layering, we obtain a 

second approximation to the actual DZ curve by utilizing the differences between the ordinates of the calculated (VES') and reference 

VES curves. As shown in step 4, of figure 4, for points on the calculated 

VES' curve that lie above the reference VES curve, we calculate new 

points (which are depicted by the arrow heads in step 4) which lie below 

the reference VES curve, and conversely. In step 5, the layering obtained from the inversion of this second approximation to the DZ curve 

is compared with the layering of the actual three -layer model. The 

calculation of the corresponding VES curve indicated that the two 

geoelectric sections are equivalent in the VES domain (that is, the two 

VES curves practically coincided). Had it been necessary to use several 

successive approximations to the DZ curve, the following iterative 

formula

_ fr° Pmi (15) 

Pmi+i Psi

would have been used where p m . +1=the (i+1) DZ resistivity, 

ps =the ordinates of the reference VES curve, p m .= ordinates of the 

DZ resistivity for the ith iteration, and the ps .=ordinate of the ith 

calculated VES curve. 

To summarize, for the first approximation we use

Pm = Pso, (16)

for the second approximation, the iterative equation 15 reduces to

and for the (i+1) to the nth approximation we use equation 15.

The simplicity of the above procedure makes it very attractive for the 

automatic interpretation of VES curves; nevertheless, the procedure 

has two drawbacks. First, it can be applied only to VES curves on which 

there are no negative slopes that are equal to or less than 1, and second, it cannot be applied to VES curves that have positive slopes that 

are equal to or slightly greater than +1. There are several types of VES 

curves which have steeply descending branches with slopes of less than 

 1, and there are also field VES curves whose slopes may be slightly 

greater than +1 because of measurement errors or lateral in-
INTERPRETATION OF SCHLUMBERGER SOUNDING CURVES E13

homogeneities. Furthermore, if the coordinates of two successive points 

on a branch of a VES curve, whose slope satisfies the relation 1 > 

slope > +1, are used in equations 13 and 14, infinite or imaginary 

resistivities and thicknesses will be obtained. In order to overcome these 

difficulties we use MDZ (modified Dar Zarrouk) equations.

MODIFIED DAR ZARROUK CURVES OF THE L TYPE

Figure 5 shows a theoretical three-layer VES curve of the H-type 

(pi>p2<ps), its corresponding DZ curve, and the Hummel curve 

(Zohdy, 1965) for the conductive second layer. The steeply descending 

branch of the VES curve occupies an intermediate position between the 

descending DZ and Hummel curves. Therefore, for steeply descending 

branches on VES curves, there are curves, here called modified Dar 

Zarrouk curves of the L type (MDZ-L), which can be calculated from 

the combination of the DZ equations and the Hummel equations and 

which will approximate that steeply descending VES branch and lie 

close to it. Furthermore, when these MDZ-L curves are solved for layering, they will result in a geoelectric section which is equivalent to the 

section for which the VES curve was calculated or observed.

The parametric equations for MDZ-L curves are

PmL =PL

'1 ^ \PLj

and

(18)

^x

where

(19)

j = average longitudinal resistivity, or Hummel resistivity, for 

the layers from 1 to j,

LL . = 2/ii= summation of layer thicknesses, or Hummel depth, to

1 the bottom of the y'th layer, 

X = number that lies in the interval 0<X<1, 

Si = hi/ PI = longitudinal unit conductance of the ith layer, 

TI =hi pi= transverse unit resistance of the ith layer. 

For X=l, equations 18 and 19 reduce to the DZ equations 11 and 12, 

and for X=0 they reduce to the Hummel equations. For values of 

0<X<1, the MDZ-L curves occupy intermediate positions between the 

corresponding Hummel and DZ curves. The upper section of figure 6
E14 TECHNIQUES IN DIRECT-CURRENT RESISTIVITY EXPLORATION

400

>-

K 100

H CC

CO hcn u<

UJ 2

CC '

H H100

52 >

en -

LU H

OC CO

o o

CC Z

uj O 10

CO _l

I III I I I I I I T I I I I

100

Hummel-\

I i i i i i

2 10 100 1000 

SCHLUMBERGER ELECTRODE SPACING (AB/2). DZ DEPTH(L m ) 

HUMMEL DEPTH (L^,). TRUE DEPTH (D), IN METRES

FIGURE 5. Three-layer VES curve and corresponding DZ curve. Hummel curve is 

drawn for the conductive second layer only. Numbers in bar designate true resistivities 

in ohm-metres.

shows a set of MDZ-L curves computed for p2/pi=0.01 and various 

values of X. The lower section of figure 6 shows a set of MDZ-L curves 

calculated for X=0.5 and various values of the resistivity ratio pt/pi.

INVERSION OF MDZ-L CURVES

The analytical expressions required for calculating the thickness and 

resistivity of a layer represented by two points on an MDZ-L curve are
INTERPRETATION OF SCHLUMBERGER SOUNDING CURVES £15

0.01

100 1000

0.1 0.5 1.5 

0.2 1 23457 10 1520 50 100 = /J 2 //I,

L-L I I ' I 'lj I 'I ' ' 1 ' ' "

1000

FIGURE 6. Modified DZ curves of the L type. Top graph: MDZ-L curves for constant 

resistivity ratio pa/pi and variable exponent X. Bottom graph: MDZ-L curves for 

constant value of the exponent X=0.5 attd variable resistivity ratio pjp\. 

hi/hi, thickness ratio; pmL, MDZ-L resistivity; LmL, MDZ-L depth.
E16 TECHNIQUES IN DIRECT-CURRENT RESISTIVITY EXPLORATION 

derived as follows. Dividing equation 19 by equation 18 we get:

mLj

which indicates that the value of S is preserved on MDZ-L curves (as it 

is preserved on Hummel and on DZ curves), and multiplying equation 

18 by equation 19 we get

h

X

(21)

For 7 = 1 in equations 20 and 21, it can be shown that

=hi, (22)

indicating that the coordinates (LmL , p mL ) of the first point represent 

the thickness and the resistivity of the first layer, respectively. In a 

multilayer section, this first layer may be a fictitious layer that replaces 

all the overlying layers. For 7 =2 in equation 20, it follows, from equation 

22, that

(23)

The thickness of the second layer, hz, can be calculated from equation 

23 if its resistivity, p2 , has already been determined. The value of p2 is 

calculated as follows. Substituting equations 22 and 23 in equation 21 

for n=2, and simplifying to get

E

-2 

^^7 L

JL

(24)

where

4 and Li+ i=L mi and L m . = abscissas of two successive H L j+i

/3y and PJ+I = PmL .and PmL . .

;

points on an MDZ-L curve, 

=ordinates of two successive 

points on an MDZ-L curve,

j+1 j+ -j, 

P2=true resistivity of the second layer.

Equation 24 contains two unknowns: X and p%. In order to solve for pz 

at a given value of 0<X<1 we use the iterative method known as the 

regula falsi method (Grove, 1966; Hildebrand, 1956).
INTERPRETATION OF SCHLUMBERGER SOUNDING CURVES E17

The successful application of the regula falsi method relies on the 

ability of the user to assume two values for p 2 , which, when consecutively substituted in the left-hand side of equation 24 here 

denoted by/ produce two values for/ with opposite algebraic signs. A 

negative value for / indicates an underestimation of p2, and a positive 

value for / indicates an overestimation. Inasmuch as the resistivities in 

a physically acceptable model must be positive, the value of zero for pi 

must be made to represent an underestimation of the true value of p2. 

Therefore, with X=l (DZ equation) and p2=0 in equation 24, if the 

value of / is positive, then the value of / must be forced to become 

negative by decrementing the value of the exponent X, from 1 to 0. 

With X=0 and P2=0, the value of /will always be negative as it can be 

readily proved from equation 24. As shown by the comparison of figures 

5 and 6, MDZ-L curves for X=0.5 approximate steeply descending 

VES branches better than DZ curves (X=l) or Hummel curves (X=0). 

Therefore, to approximate descending VES branches, the "first" value 

of X is generally taken as 0.6, and, if /is positive, then X is decremented 

by 0.2 (X=0.6, 0.4, 0.2, 0). The first value of X and the rate at which it 

is decremented affect the thicknesses and resistivities of the layers to be 

obtained from the inversion.

By assuming p 2 =Pi+i, a positive value of/, /+, always is obtained 

because, on a descending MDZ-L branch (for any value of 0<X<1), 

the MDZ-L resistivity pi+l always is greater than p2 except at the 

asymptote where fc+i=p2, if #+1=02, then the value of / is zero, 

and the required value of pz is determined.

Having determined the first value for X=X' (where 0<X'<0.6), for 

which /=/_ at P2=P2-=0, and using p 2+=p'j+ i to calculate /=/+ at 

X=X', we use the iterative regula falsi equation to calculate the first 

approximation for pz from

(25)

We substitute pz (1) for p2 in equation 24 to calculate /_ (1), then we 

substitute p2 (1) and /_ (1) for p 2 and f_ in equation 25 to calculate the 

second approximation for pz=pz (2). The iteration is repeated for a 

maximum of 15 times or until

0.98<_£2(iL<l.02. (26) P2 fl-D

Having determined the value of p2 , by the above iterative method we 

can calculate the thickness, hz, from equation 23.

To summarize, for the automatic interpretation of VES curves, 

ascending branches on VES curves are approximated by DZ curves 

and descending branches are approximated by MDZ-L curves. The
E18 TECHNIQUES IN DIRECT-CURRENT RESISTIVITY EXPLORATION

iterative formula 15 is used for calculating the ordinates of the successive DZ and MDZ-L points whose inversions successively approximate 

the required solution.

MODIFIED DAR ZARROUK CURVES OF THE T TYPE

A theoretical VES curve for horizontally stratified, laterally 

homogeneous media may have a rising branch with a slope of +1, which 

reflects the detection of a layer with an infinitely large resistivity; in 

practice, field VES curves may have slopes that are slightly, or 

significantly, greater than +1. Such slopes may be caused by lateral 

heterogeneities (Kalenov, 1957), fences with metal posts, errors in 

measurements, or current leakage (Zohdy, 1968). Furthermore, in 

applying the iterative equation 15 for adjusting the ordinates of DZ and 

MDZ-L curves, it is possible that the slope between two successive 

points may become >1 after several iterations have been made. If the 

ordinates of points on such steeply ascending segments are used in the 

DZ inversion equations, infinite or imaginary resistivities and 

thicknesses will be obtained.

To overcome this problem, we define MDZ-T curves whose 

parametric equations are given by

Pmtj- M- "~Pt'W) '

and

Sh,

(27)

where

(28)

pt .=27yz/i,=average transverse resistivity of the layers from 1 

to;,

Lt .= 2/it= summation of layer thickness to the bottom of the y'th 

layer,

X=number that lies in the interval 0<X<1.

Equations 27 and 28 are similar to equations 18 and 19 except that ptt 

instead of PL, is used to modify the DZ resistivities and depths. 

Depending on the value of X and on the parameters of the section, the 

corresponding MDZ-T curve can have slopes that are greater than 1. In 

fact, the MDZ-T curves (fig. 7) look almost like a mirror image of 

MDZ-L curves (fig. 6) reflected across the abscissa.
INTERPRETATION OF SCHLUMBERGER SOUNDING CURVES E19

100

1000

100

s 1 °-

iiii.ii i .1 . . i..

100 1000

FIGURE 1. Modified DZ curves of the T type (MDZ-T). Top graph: MDZ-T curves for 

constant resistivity ratio p2/p\ and variable exponent X. Bottom graph: MDZ-T 

curves for constant value of the exponent X=0.5 and variable resistivity ratio PZ/PI. 

hi/hi, thickness ratio; p mt, MDZ-T resistivity; L mt , MDZ-T depth.
E20 TECHNIQUES IN DIRECT-CURRENT RESISTIVITY EXPLORATION

INVERSION OF MDZ-T CURVES

The equations for the inversion of MDZ-T curves are similar to those 

for the inversion of MDZ-L curves and are derived as follows. Multiplying equation 27 by 28 we get

(29)

which indicates that the values of T are preserved on MDZ-T curves. 

Dividing equation 27 by equation 28 we get

j i j_ i *^

P'"" ~" """ (30)

For i = l in equation 27 and 28, it can be shown that

pm =piandL =/iu (31)

and for 7=2 in equation 30, it follows, from equations 29 and 31 that

P2

(32)

substituting equation 31 and 32 in equation 30 for j=2, and simplifying 

we get

ft

-1=0, (33)

where

= L.andL mt m

and -= and

P2=true resistivity of the second layer.

Equation 33 is analogous to equation 24 and can be solved for p2 , using 

the regula falsi method. For an underestimation of pz we set pz=Pi+i, 

and for an overestimation of p2, we set p2=50/5J+i. The equation is 

tested with p2= 50 p^fi and X=l, 0.9, 0.8, . . . and so on, in search of 

a positive value for the left-hand side, which would indicate that p2=50 

#+! is an overestimation. Then, having guaranteed that the appropriate value of X has been determined for a positive value of p2 , we
INTERPRETATION OF SCHLUMBERGER SOUNDING CURVES E21

proceed to use the regula falsi method for evaluating p2 and, hence, hz 

from equation 32.

INVERSION OF VES CURVES

The inversion of a Schlumberger VES curve is made by finding an 

MDZ curve (which may be formed of a combination of DZ, MDZ-L, 

and MDZ-T segments) which lies close to the observed VES curve and 

whose inversion results in a layering for which a calculated VES curve 

practically coincides with the observed VES curve. The slopes at the 

successive points on the digitized VES curve are calculated using the 

formula

SLOPEi+1 =\og JB.. /»../ log Ll+1 Lt , (34)

where L= AS/2 = Schlumberger electrode spacing, and i=l, ... , n. 

Equation 34 assigns the value of the slope between the ith and (i+1) 

point to the (i+1) point. The slope of the first point is assumed to be 

equal to zero. 

There are at least two advantages for the computation of these slopes:

1. If the slope exceeds +1, the computer will print a warning message 

about the distortion of the VES curve, and if the slope exceeds the arbitrarily set value of 1.4, then this generally indicates that an error in 

punching or digitizing the data has been made, and the problem is 

automatically rejected by the computer unless it was specifically instructed to find the best fitting curve regardless of curve distortions.

2. The fitting tolerance, FT, for each point on the observed VES 

curve, is defined in terms of the slope by the formula

FTi=M+N (SLOPEi) 2 , (35)

where M and N are usually taken equal to 5 and 1, respectively. The fitting tolerance for the first point is arbitrarily assigned the value of 5. 

Equation 35 allows a larger fitting tolerance for steeply descending 

branches than for moderately descending or ascending branches.

The computer compares the values of the successive ordinates on the 

digitized VES curve and, for each pair of successive points, it 

decides, on the basis of the value of Pm,- +1 /Pmj> whether the DZ, 

MDZ-L, or MDZ-T inversion equations should be used. By applying 

the appropriate inversion formulas, we obtain a geoelectric layering in 

which the number of layers is equal to the number of points on the 

digitized VES curve. The thickness of the last layer is set to equal a very 

large number (9999999).

The TKF curve is calculated for this layering using Sunde's recursion 

formula and a VES curve is calculated by convolution using Ghosh's 

inverse filter coefficients, so that the abscissas of the calculated and observed VES points are the same. The percentage
E22 TECHNIQUES IN DIRECT-CURRENT RESISTIVITY EXPLORATION

difference, PD, between calculated pc , and observed, po, apparent 

resistivity values is calculated from

X100

1 log p0 . 

The sum of the squared residuals (SSQR)

SSQR= 2

i=i

is also calculated.

If the percentage difference for one calculated point exceeds its 

precalculated fitting tolerance, the iterative equation 15 is used to 

calculate the ordinates of a new set of MDZ points. This newly constructed MDZ curve is inverted, and a second VES curve is calculated. 

If a match between observed and calculated VES curves (within the 

prescribed fitting tolerance for each point) is not obtained after 10 

iterations (using equation 15), or if the sum of the squared residuals increases 5 times, instead of decreasing, then the layering solution with 

the least SSQR is saved in the computer.

The failure of the iterative procedure to produce a satisfactory match 

between observed and calculated VES curves indicates one or both of 

the following possibilities:

1. The curve is distorted and therefore it may be impossible to fit every 

point on it with a theoretical curve, which is calculated for 

horizontal layering.

2. The left branch of the curve is incomplete, and (because, in the 

iterative procedure, only the ordinate of the first point is allowed 

to change and not its abscissa) it is impossible to fit the curve with 

a model in which the thickness of the first layer is equal to the 

electrode spacing of the first point. This problem can be solved if 

the thickness of the first layer is forced to become less than the 

abscissa of the first point on the curve.

Therefore, we introduce the concept of the thickness reduction factor 

0<y<l, which is multiplied by the value of the thicknesses h\ and hz 

after their calculation from equations 14, 23, or 32. Thus, for the first set 

of iterations, the value of the thickness reduction factor, y, is set equal 

to unity. If this does not result in a satisfactory match, then, with the 

value of y decreased by an decrement of 0.1, the whole iterative 

procedure is repeated from the start. For each set of iterations, with a 

fixed value of y, the least SSQR obtained is compared with the one obtained from the previous set of iterations. If the value of the least SSQR 

for one set of iterations becomes larger than that of the previous set of 

iterations, the automatic inversion procedure is terminated and the 

solution with the least SSQR is printed. The maximum number of
INTERPRETATION OF SCHLUMBERGER SOUNDING CURVES E23

allowed iterations is set to equal 60 (which implies that the value of y is 

allowed to decrease a maximum of 5 times).

Figure 8 shows the output of five iterations from the first set of 

iterations (in which the value of y is equal to unity) for the automatic 

interpretation of a theoretical three-layer curve of the H type. In step 1 

the sought MDZ curve is approximated by the observed VESo curve 

itself (MDZi = VESo). The inversion of this MDZi curve and subsequent calculation of the apparent resistivities result in the VESi curve. 

The ordinates of the VES0 and VESi curves are compared and the ordinates of the second approximation (MDZ2 ) for the MDZ curve are 

calculated from equation 15 which is symbolized differently (from 

equation 15) on the figure (pSQ = VES0 , ps ., = VES,-, pm . = MDZ,-, and 

Pmi+ l = MDZ,-+1). The iteration is continued through step 5 where 

the curve of VESs, which is calculated for the layering obtained from 

the inversion of the MDZs curve, practically coincides with the curve of 

VESo. In this example, VESi through VESs were calculated by 

numerical integration rather than by convolution.

The number of layers in the resulting model is equal to the number of 

points on the digitized curve. This model constitutes the detailed solution, which, for the example under consideration, is composed of 15 

layers. In the detailed solution, some successive layers may have almost 

the same resistivity and, therefore, the actual number of layers with 

distinctly different resistivities may be less than the number of points 

on the curve. In this example the number of such distinct layers is nine. 

It is important to realize that had the reference VES curve been 

calculated or observed over these nine layers, then its interpretation, by 

conventional curve-matching methods, would have been made in terms 

of the three-layer model.

Figure 9 shows the DZ curves and the layering distributions for both 

the automatically obtained nine-layer model and the original 

theoretical three-layer model. The resistivities of the infinitely thick 

bottom layers on the two models are different because the last branch 

on the VESo curve is not well developed and was not plotted to the point 

at which it asymptotically would have approached the value of the true 

resistivity of the bottom layer. Nevertheless, this example illustrates 

the fact that the application of the method does not depend on the extension and full development of the last branch on the observed VES 

curve, which is a requirement for the other methods of automatic interpretation that use the transformation of VES curves into their corresponding total-kernel-function curves by convolution (Strakhov and 

Karelina, 1969; Ghosh, 1971a). The development of the first branch on 

the observed VES curve, at least to electrode spacings that are equal to 

or less than the thickness of the first layer, is recommended but not required.
E24 TECHNIQUES IN DIRECT-CURRENT RESISTIVITY EXPLORATION

100 

100

100

100

100

10

Layerin9 MDZ- = VES0 X MDZ, 

VES.

©

I I I I 111 I I I I I I

10 100 1000 

AB/2,L MDZ ,\N METRES

10,000
INTERPRETATION OF SCHLUMBERGER SOUNDING CURVES E25

<1 FIGURE 8. Curves showing the output of five iterations for the inversion of a theoretical 

three-layer VES curve of the H type. VES0 , reference VES curve or ordinates on VES0 

curve; MDZi to MDZs, successive approximations to ordinates of MDZ curve; VESi to 

VESs, VES curves calculated for the layering obtained from the inversion of 

MDZi to MDZ5 curves, respectively. p MDZ , MDZ resistivity; L MDZ, MDZ depth; ps, 

Schlumberger apparent resistivity; AB/2, Schlumberger electrode spacing.

400

100

100

0.4

3-leyer model 

9-layer model

i i i i i i i i i i i i i

10 100 1000

DZ DEPTH (Lm ), AND DEPTH (D), IN METRES

10,000

FIGURE 9. Comparison between equivalent DZ curves and equivalent layering for 

reference theoretical three-layer model and automatically obtained nine-layer model.
E26 TECHNIQUES IN DIRECT-CURRENT RESISTIVITY EXPLORATION

INVERSION OF INCOMPLETE AND (OR) 

DISTORTED VES CURVES

Incomplete VES curves are here defined as VES curves on which 

neither the left branch nor the right branch forms asymptotes to the 

first- and last-layer resistivities, respectively. However, a VES curve 

whose right branch rises at an angle of 45° and is defined by at least 

three points extending over about one-third of a logarithmic cycle is 

considered as having a complete terminal branch.

In contrast to other automatic or semiautomatic methods of interpretation, the method presented here does not rely on the completeness of the first and terminal branches of the observed curve. The 

problem created by the incompleteness of the left branch is overcome 

by introducing the thickness reduction factor, y, mentioned earlier, and 

the problem of the incompleteness of the right VES branch is simple to 

resolve by using the MDZ method and setting the thickness of the last 

layer equal to a very large number. However, when the right branch on 

the VES curve is incomplete, the automatically interpreted model may 

show a resistivity for the last layer which is different from the resistivity 

anticipated by the interpreter; nevertheless, all the given points on the 

observed curve would be properly fitted. Therefore, if the interpreter 

has any preconceived notions about the resistivity of the last layer, then 

he should extrapolate the terminal branch toward its assumed asymptote.

Although the completeness of the first (left) branch is not a requirement for the applicability of the computer program (Zohdy, 1974c), it 

simplifies the problem, and a few seconds of computer time may be 

saved.

Distorted VES curves are defined here as VES curves whose curvature at certain points cannot be fitted with theoretical VES curves 

which are calculated for horizontally stratified media. The distortion on 

a VES curve may be caused by inhomogeneities in the horizontal direction or by errors in measurement or may even be created by errors in 

digitizing the data for an otherwise smooth VES curve. If the distortion 

on a given VES curve is such that the logarithmic slope on the curve at 

any given point exceeds 1.4, then, as mentioned earlier, the computer 

automatically rejects the problem unless it was instructed not to. If, 

however, the distortion on the VES curve is more moderate, then the 

automatic interpretation program will yield a calculated VES curve 

which fits the digitized points in the least-square sense, although not 

necessarily within the prescribed fitting tolerance.

In attempting, even though unsuccessfully, to fit all the given points 

on a distorted VES curve, the detailed layering may contain a 

geologically improbable sequence of thin layers with highly contrasting 

resistivities. Consequently, the calculated VES curve for the
INTERPRETATION OF SCHLUMBERGER SOUNDING CURVES E27

improbable layering (which fits the distorted VES curve in the leastsquare sense) is considered as a smoothed version of the observed curve 

and is automatically reinterpreted to obtain a more realistic layering.

In the reinterpreted model, the replacement of the sequence of layers 

having highly contrasting resistivities by a sequence of layers having 

more moderately contrasting resistivities shows that the layers with 

highly contrasting resistivities were not only geologically improbable 

but also geoelectrically unnecessary; if their presence was accurately 

reflected on the calculated VES curve before reinterpretation, then, 

most likely, they would have been obtained also in the reinterpreted 

model.

The computer program is written so that if one point on the digitized 

observed VES curve cannot be fitted within the prescribed fitting 

tolerance, then the detailed solution for which 2[logpo logpj 2 = 

minimum is chosen. The VES curve for this detailed model is 

recalculated, starting at an electrode spacing equal to one-tenth the 

smallest given electrode spacing and, if specified by the user, ending at 

an electrode spacing equal to ten times the largest given electrode spacing (thus generating a complete VES curve if the given curve was incomplete as well as distorted). This new VES curve, which is complete 

and smooth, is automatically reinterpreted. The number of layers in the 

detailed model is reduced through the automatic smoothing of the corresponding DZ curve.

Figure 10 shows an example of the interpretation of an incomplete 

and distorted VES curve. The layering shown in figure 10A contains an 

unrealistic sequence of thin layers, with resistivities ranging from 1.05 

to about 2,600 ohm-metres. These layers were created primarily in attempting to fit the incomplete left branch of the curve and the distorted 

segment on the VES curve between the electrode spacing of 100 and 

1,000 m (metres).

Figure 10B shows the DZ curve for the unrealistic detailed model. 

The fundamental DZ points occurring at the DZ depths of about 75, 

110, and 172 m, respectively, cannot be easily smoothed by using the 

automatic procedures to be described later, although they can be 

smoothed manually.

Figure 10C shows the complete VES curve (starting at the electrode 

spacing of AB/2 = 1 m and terminating at AB/2 = 10,000 m calculated 

for the unrealistic detailed model shown in figure WA. This VES curve 

(which now may be considered as the complete and smoothed version of 

the observed curve) was automatically interpreted in terms of the 

realistic detailed layering also shown in figure 10C; the VES curve 

calculated for this layering is shown in the graph with the open circles.

Figure IQD shows the detailed DZ curve for the 25 layers in figure 10C 

and the fundamental DZ points (1-6) which were obtained using the
E28 TECHNIQUES IN DIRECT-CURRENT RESISTIVITY EXPLORATION

CO

Ul

cc

Ul

^

7 1000

I

O

Z

§

> 100

H

co

CO

Ul

cc

UJ

D

CC

1 10

< 100

H 

H

en 10 

uj 100

CC

N

O

"j

>

1 li8

CO

Ul

cc

H

z

UJ

cc<0.

< 10 

< 100

cc

Ul

o

cc

Ul

m

D

_J

I

O 

co m

i i i i i i i 1 1 i i i i i i 1 1 i i i i i i 1 1 1 i i i i i i i i

-

-

-

_

-

-

_

A

"

i i

 -

-

B /

-

i 

C

13--°'

rn

| |

*-,

\

"

-

-

-Detailed model for distorted VES curve

^sZ* ~

/S^ Observed curve

^>*» -9^^ Best fitting VES curve

~

-

 "^DZ curve for detailed model in A

1 1 ^o ° °" ° °

^-°

o^"" °^ ^° Smoothed and _

a/

^°

^ _

°--~ f. x^'0 completed curve

 ^Calculated curve

^Detailed layering for smoothed curve

- 

: 3 JS-^-""^ '~"~~'~t

- °- »--5s^4 *^^ * 

x'" ^^^\ .^^*

_. / v /S*~~ Detailed DZ curve

 ^ 2/* ^» ^x for model in C ? >JL 

1 j* 5^6-layer equivalent DZ curve

 -*=&

-

: 2

-

#

lo _,o"

3 

6

4 ^**~°"°"<>

^ 0 ''S

) 0 O

^"*\ _^ff^

"~*'^^:?z^6- layer VES curve

_

 6-layer model

o

1 1 1 1 M 1 1 1 II II 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

1 10 100 1000 10,000 

ELECTRODE SPACING (AB/2). DZ DEPTH (Lm ). AND DEPTH (D). IN METRES

FIGURE 10. Graphs showing the automatic interpretation of an incomplete and

distorted VES curve.
INTERPRETATION OF SCHLUMBERGER SOUNDING CURVES E29

automatic-DZ-smoothing subroutine (Zohdy, 1974c). The automatic 

inversion of the smoothed DZ curve resulted in the six-layer model 

shown in figure IQE, where the corresponding six-layer VES curve is 

compared with the original digitized data of the observed curve.

AUTOMATIC SMOOTHING AND INVERSION OF DZ CURVES

The number of layers in the detailed model obtained from the 

automatic interpretation can be reduced automatically so that the 

resulting model will contain only the fundamental layers in the geoelectric section, and so that the reduced model will be electrically 

equivalent to the detailed model. The number of layers is reduced by 

first calculating the DZ curve for the detailed model, then smoothing it 

numerically and inverting it. The following two methods were devised 

for smoothing the detailed DZ curve.

1. FIRST METHOD

In this first method the fundamental DZ points designating the termination of one fundamental DZ branch and the beginning of another 

are also points on the detailed DZ curve. By successively substituting the coordinates of the first DZ point (Lmi, pmi), which is taken 

as the first fundamental DZ point, and the coordinates of the successive 

points (Lmt , Pmt ), (Lma> Pm3 ),..., (Lmk, Pmk ) in the inversion equation

I /Lmk Pmk Lmi pmi

II Lmk _ Lmi 

* o , a _.

P2* k =f I 1 2 ^_', (36)

Pmk

we can solve for a succession of second-layer resistivities, p2*. Each time 

a value for p2 * is calculated, we use the equation (Zohdy, 1974a)

(37)

to calculate p m values at the abscissas L m . which lie between L mi 

and Lmk. We calculate the ratio pmc /pm to evaluate the deviations of 

the ordinates on the calculated two-layer DZ curve from the ordinates (p m) of points on the detailed DZ curve. If, upon using the coordinates of the (/+!) point on the detailed DZ curve, the 

ratio Pmc /Pm, for any point between 1 and (7+1), exceeds 1.05 or 

becomes less than 0.95, the (/+!) point is dropped and the coordinates 

of the ith point on the detailed DZ curve are taken as the coordinates of 

the second fundamental DZ point.
E30 TECHNIQUES IN DIRECT-CURRENT RESISTIVITY EXPLORATION

The above procedure is repeated with the last-determined fundamental DZ point taken as the starting point. When tlie nth point on the 

detailed DZ curve has been used, the smoothing procedure is terminated, and the coordinates of the fundamental DZ points are used in 

the inversion equations 13 and 14 to calculate the reduced layering.

2. SECOND METHOD

This method is based in part on applying the method of least squares 

to the linear representation of DZ curves, which is a plot of T as a function of S. For a homogeneous and isotropic earth, the curve T=f(S) is a 

straight line. This can be shown by considering the following equations. 

For an inhomogeneous earth,

T=Plhi+p2h2+. . .+Pn hn= Pt H, (38) 

S=tll.+JlL+. . .+£s=#, (39)

Pi P2 Pn pL

T=ptpLS, (40) 

but, for a homogeneous earth, p t=pi=p, and, therefore,

T=PiS, (41)

which is the equation of a straight line whose slope is equal to the 

square of the resistivity of the medium.

In bilogarithmic presentation the first fundamental DZ branch 

(pm =f(Lm ))is a horizontal line; preferably the smoothing process for 

the first fundamental branch is made by successively averaging the 

ordinates of the successive DZ points on the n-layer DZ curve, using the 

formula

Pm*(D=e , (42)

where

i = l, 2,. .... ^successive number of points,

/= number of averaged points, and

p m *(I) = aver age value of pm in the logarithmic sense. 

The successive deviations

5,(/)=Pm P.*(/) (43)

are calculated, and if with the inclusion of the (/+!) point the value of 

6,(I+1), for any value of l<i</+l, is greater than 1.05 or smaller than 

0.95, the averaging process is terminated, and the value ofpm(I) is taken 

as equal to p\, which is the average true resistivity for the "first" layer.
INTERPRETATION OF SCHLUMBERGER SOUNDING CURVES E31

In linear presentation of DZ curves, the general equation for the first 

straight line is

Ti=A Ql +A nSi, (44)

where Aoi=0, and An=pi*'2 (see equation 41). Starting at the last point 

that was included in the above logarithmic averaging process, we 

successively use the least-squares equations

2T-SS (45)

(46)

where T= 2T, , S- 2SZ- , and /= number of points used, to compute the 

successive coefficients Ao and A\ of a succession of 

the straight lines:

). (47)

The abscissas (S) and ordinates (T) of the intersections of these 

successive straight lines with the first straight line are calculated 

successively from the equations

Si = A02 (7)-Aoi (48) 

An-Ai 2 (D '

T^Aoi+AnSi. (49)

Equations 48 and 49 are obtained by the simultaneous solution of 

equations 44 and 47 for Ti = T2 and S\ = S2 .

The number of points, /, to be included on any given T=f(S) straight 

line or DZ branch is governed by the following conditions:

1. If, upon the inclusion of (/+!) point for the Jth line, the value of the 

slope (A\j) becomes negative or zero, an imaginary or zero 

resistivity (pf = ^A\j) is obtained and the (/+!) point is dropped.

2. If, upon the inclusion of the (/+!) point on the Jth line, the intersection of the Jth and (J 1) lines occurs at an abscissa value less 

than the abscissa of the point preceding the last point on the (J 1) 

line or greater than the abscissa of the point succeeding the last 

point on the (J-l) line, then the (/+!) point is dropped.

3. If, in condition 2 above 7+1 = 2, then the coordinates T and S of the 

first point on the Jth line (which are also the coordinates of the last 

point that was averaged on the J 1 line) are redefined in terms of 

T and S of the corresponding point on the smoothed (J 1) line. 

This will guarantee that the coordinates of the point of intersection of (J-l) line and the Jth line are the same as the coordinates 

of the last point on the smoothed (J-l) line.
E32 TECHNIQUES IN DIRECT-CURRENT RESISTIVITY EXPLORATION

4. If, upon the inclusion of the (/+ 1) point, the fitting tolerance for any 

point 1 <i<(/+l) is exceeded, then, again, the (/+!) point is 

dropped.

The fitting tolerance (in condition 4) is here defined in terms of the 

slopes, si, of the points on the two-layer DZ curves (bilogarithmic 

presentation) which are used to fit the fundamental segments on the nlayer DZ curve. The following equation is used to define the fitting 

tolerance:

, (5Q)

where

a = a number set equal to 2, 4, or 8, 

s,=slope of the ith point on the two layer DZ curve. 

A plot of the function K",=/(SJ) for a = 1.5, 2, 2.5, and 3 is shown in figure 

11 and indicates that the fitting tolerance generally is larger for positive 

values of s,- than for negative values of s( . This property makes the 

function K suitable for defining the fitting tolerance in accordance with 

the rules governing the equivalence between DZ curves and their corresponding VES curves (Zohdy, 1974a). Furthermore, the value of K 

is always greater than or equal to 1 for all values of s,- in the range 

 l<Sj<l, which makes it suitable for comparison with the values of 

the deviations

ti(D=Pm/Pmc - (51)

Thus, for each set of points on the T=f(S) line which are fitted by the 

method of least squares, an intercept with the previous straight line is 

determined, and the values of T and S at the intersections are used to 

define the coordinates of pm *=p\*, L^=h\*, from the equations

L

and p2* is calculated from the slope of the second straight line

T*=f(S*).

The values of pm are calculated successively from equations 37 at the

abscissas, Lm, for all the averaged DZ points, and if upon the inclusion 

of the (/+!) point, the deviation <5,(/) for l<i</+l is greater than 

the calculated value of K( or smaller than 1/K{ , the (/+!) point is 

dropped.

The above procedure is successively applied to evaluate the coordinates (Lm *, pm*) of all the fundamental DZ points on the n-layer DZ 

curve. Then, the inversion equations 13 and 14 are used to determine 

the layer thicknesses and resistivities of the reduced model. The VES
INTERPRETATION OF SCHLUMBERGER SOUNDING CURVES E33

-1 -0.5 -0.25 0 0.25 0.5

DZ SLOPE

FIGURE 11. Graph showing the variation of the fitting tolerance, K, as a function of the 

slope of a point on a DZ curve. Values of a in equation 15 are shown on each curve.

curve for the reduced model is calculated to ascertain the equivalence 

between the calculated and observed models. In the computer program 

the maximum number of layers in the reduced model is set equal to 10. 

Thus, if upon the completion of the smoothing process the resulting 

number of layers is larger than 10, then the constant a in the fittingtolerance function K is increased from 2 to 4, and, if necessary, to 8, and 

the smoothing process is repeated.

FIELD EXAMPLE

Figure 12 shows an example of the results obtained from the 

automatic interpretation of a deep vertical electrical sounding, VES 19, 

obtained by Zohdy and Stanley (1973) on the Snake River plain, Idaho. 

The observed VES curve at the top of the figure was measured with the 

Schlumberger array to electrode spacings of AB/2 = 12,000 feet («3,657 

m) and with the equatorial array to a maximum spacing of about 68,000 

feet («20,726 m). The various segments on the curve were shifted 

downward to conform to the equatorial segment (adjusted curve), and 

the minor cusps, caused by the crossing of the current electrodes over 

lateral heterogeneities, were removed to form the smoothed curve.

As shown in the lower part of the figure, the smoothed curve was
E34 TECHNIQUES IN DIRECT-CURRENT RESISTIVITY EXPLORATION

4000

1 10 100 1000 10,000

SCHLUMBERGER SPACING (AB/2), EQUATORIAL SPACING (R),

INTERPRETED DEPTH (D). IN METRES

FIGURE 12. Results of the automatic interpretation of a deep vertical electrical 

sounding (VES 19) curve obtained on the Snake River Plain, Idaho.

digitized at a logarithmically equal interval, A, which corresponds to 

the rate of 6 points per logarithmic cycle (A = e ln 10/6 ) to facilitate the 

use of the Ghosh convolution coefficients (Ghosh, 1971b). The data 

(which consist of the smallest and largest electrode spacings (AB/2) and 

the values of apparent resistivities at the prescribed digitizing interval) 

were fed into the computer which performed the automatic interpretation and fitted the digitized curve with the calculated VES curve, as 

shown in the lower part of the figure. The detailed solution (27 layers 

corresponding to 27 points on the digitized curve) is also depicted on the 

lower part of the figure.

There are basically five electrical units on the detailed model, as 

evidenced by the five major breaks in the distribution of resistivity with
INTERPRETATION OF SCHLUMBERGER SOUNDING CURVES E35

depth. The first unit, with an average resistivity of about 30 ohmmetres, corresponds to a windblown-soil layer. The second unit, with an 

average resistivity of about 4,000 ohm-metres, represents dry basalts of 

the Snake River Group. The third unit, at a depth of about 

300-400 feet (91.44-121.9 m), is characterized by an average resistivity 

of about 300-400 ohm-metres and is interpreted as basalt saturated 

with fresh water. The fourth layer is characterized by low resistivities 

which range from about 10 ohm-metres to about 100 ohm-metres; most 

likely these resistivities represent sedimentary rocks, which are possibly 

intercalated with ash-flow tuffs, and which extend to a depth of about 

20,000 feet (6,096 m). The fifth electrical unit is the electric basement, 

which has a high resistivity (> 1,000 ohm-metres) and may represent 

Paleozoic rocks of very large thicknesses.

The number of layers in the detailed solution can be reduced either 

automatically (as explained earlier) or manually by plotting the DZ 

curve for the 27 layers, smoothing it, and inverting it. Figure 13 shows 

the 27-layer DZ curve as the dashed curve. A five-layer solution (equations 13 and 14) was obtained by manually selecting the coordinates of 

the DZ points (1, 2, 5, and 6) designated by the open circles and a sixlayer solution was obtained by manually selecting the coordinates of the

5000

330 780 44 oo Q M

792 1.400 7,300 METRES

10 100 1000 10,000 100,000 

DZ DEPTH (Lm ), INTERPRETED DEPTH {£», IN METRES

FIGURE 13. Graph showing two equivalent-layering models obtained by manually 

smoothing the DZ curve of the detailed solution (dashed curve) and inverting the 

smoothed DZ curves. Model I is based on the inversion of the DZ points 1, 2, 5, and 6. 

Model II is based on the inversion of the DZ points 1, 2, 3, 4, and 7-
E36 TECHNIQUES IN DIRECT-CURRENT RESISTIVITY EXPLORATION

first two open-circle points (points 1 and 2) and the following three 

solid-circle points (points 3, 4, and 7). Other equivalent models can be 

found by automatically or manually smoothing the DZ curve in different ways.

ADVANTAGES AND LIMITATIONS OF THE METHOD

The advantages of the method of automatic interpretation presented 

in this paper are:

1. The average processing time for a VES curve covering three 

logarithmic cycles is about 8 seconds on an IBM 360/65 computer. 

This makes the method very suitable for the processing of large 

numbers of VES curves at considerable savings in time and expenditure.

2. Only the values of the smallest and largest electrode spacings, and 

the ordinates of the apparent resistivities (which must have been 

digitized at the logarithmic interval of 1.46 along the abscissa) are 

entered into the computer.

3. The method does not rely on the extension of the left and right 

branches of the VES curve to their appropriate asymptotes.

4. Distorted field VES curves can be processed and two sets of solutions 

will result: one based on the inversion of the distorted curve and 

the other based on the inversion of the automatically smoothed 

VES curve.

5. The method is extended through the automatic smoothing of DZ 

curves so that earth models, composed of 10 layers or less, are 

automatically obtainable.

6. The user does not have to make any preliminary guesses at the 

number of layers, their thicknesses, or resistivities.

7. The curve-fitting process is made in the VES curve domain rather 

than in the kernel-function domain. This assures the interpreter of 

obtaining a valid solution which fits the observed data and is 

within the possible equivalence for VES curves rather than the 

wider range of equivalence for TKF curves. 

The limitations of the method primarily lie in the limited accuracy of

computing VES curves by convolution, using Ghosh's inverse filter

coefficients. The calculation of a large variety of VES curves, using

Ghosh's nine coefficients indicated the following:

1. The position of the S line will be shifted to the left by about 6 percent 

if the resistivity of the last layer is assumed to be much greater 

than about 100 times the overlying layer. This is not a serious 

problem because, by assuming the resistivity of the last layer to be 

about 50-100 times the resistivity of the overlying layer, a sufficiently long linear segment on the VES curve will be calculated, 

and it will coincide with the true S line within about 1 percent.
INTERPRETATION OF SCHLUMBERGER SOUNDING CURVES E37

2. For VES curves with steeply descending branches and for which the 

ratio of the minimum apparent resistivity, PSmin> to the maximum apparent resistivity, ps , is less than 0.025, the application XTlaX

of Ghosh's coefficients will result in inaccurate apparent resistivity values for values of p s . which are less than 0.025 ps "mm *max

Consequently, the automatic fitting of such curves may be less 

than satisfactory in the region where p s . <0.025 ps . For these rnin rnflx

curves the validity of the solution should be checked by calculating the VES curve for the obtained model using numerical integration methods.

REFERENCES CITED

Bhattacharya, P. K., and Patra, H. P., 1968, Direct-current geoelectric sounding Principles and interpretation: New York, Elsevier Publishers, 139 p. 

Compagnie Generate de Geophysique, 1963, Master curves for electrical sounding [2d rev.

ed.]: The Hague, European Assoc. Explor. Geophysicists. 

Crous, C. M., 1971, Computer-assisted interpretation of electrical soundings: Colorado

School Mines M.S. Thesis, 108 p. 

Flathe, H., 1955, A practical method of calculating geoelectrical model graphs for

horizontally stratified media: Geophys. Prosp. [Netherlands], v. 3, no. 3, p. 268-294. 

___1963, Five-layer master curves for the hydrogeological interpretation of geoelectric

resistivity measurements above a two-storey aquifer: Geophys. Prosp. [Netherlands], v.

11, no. 4, p. 471-508. 

Ghosh, D. P., 1971a, The application of linear filter theory to the direct interpretation of

geoelectrical resistivity sounding measurements: Geophys. Prosp. [Netherlands], v. 19,

no. 2, p. 192-217.

___1971b, Inverse filter coefficients for the computation of apparent resistivity standard curves for a horizontally stratified earth: Geophys. Prosp. [Netherlands], v. 19, no.

4, p. 769-775. 

Grove, W. E., 1966, Brief numerical methods, in Applied mathematics series: New Jersey,

Prentice-Hall, 117 p. 

Hildebrand, F. B., 1956, Introduction to numerical analysis: New York, McGraw-Hill

Book Co., 511 p. 

Kalenov, E. N., 1957, [Interpretation of vertical electrical sounding curves]: Moscow,

Gostoptekhizdat, 471 p. [in Russian]. 

Koefoed, Otto, 1965, Direct methods of interpreting resistivity observations: Geophys.

Prosp. [Netherlands], v. 13, no. 4, p. 568-591. 

 1966, The direct interpretation of resistivity observations made with a Wenner

electrode configuration: Geophys. Prosp. [Netherlands], v. 14, no. 1, p. 71-79.

.1968, The application of the kernel function in interpreting geoelectrical measurements: Geoexploration Mon., ser. 1, no. 2, Stuttgart, Gebriider Borntraeger, 111 p. 

-1970, A fast method for determining the layer distribution from the raised kernel

function in geoelectrical sounding: Geophys. Prosp. [Netherlands], v. 18, no. 4, p. 

564-570.

Kunetz, Geza, 1966, Principles of direct-current resistivity prospecting: Berlin, Gebriider 

Borntraeger, 103 p. [English translation from the French by Robert Van Nostrand]. 

Kunetz, Geza, and Rocroi, J. P., 1970, Traitement automatique des sondages electriques: 

Geophys. Prosp. [Netherlands], v. 18, no. 2, p. 157-198.
E38 TECHNIQUES IN DIRECT-CURRENT RESISTIVITY EXPLORATION

Maillet, Raymond, 1947, The fundamental equations of electrical prospecting:

Geophysics, v. 12, no. 4, p. 529-556. 

Meinardus, H. A., 1967, The kernel function in direct-current resistivity sounding:

Colorado School Mines D. Sc. Thesis, 151 p. 

___1970, Numerical interpretation of resistivity soundings over horizontal beds:

Geophys. Prosp. [Netherlands], v. 18, no. 3, p. 415-433. 

Mooney, H. M., Orellana, Ernesto, Pickett, Harry, and Tornheim, Leonard, 1966, A

resistivity computation method for layered earth models: Geophysics, v. 31, no. 1, p.

192-203. 

Onodera, Seibe, 1960, The kernel function in the multiple-layer resistivity problem: Jour.

Geophys. Research, v. 65, p. 3787-3794. 

Orellana, Ernesto, 1963, Properties and drawing of the so-called Dar Zarrouk curves:

Geophysics, v. 28, no. 1, p. 99-110. 

Orellana, Ernesto, and Mooney, H. M., 1966, Master tables and curves for vertical

electrical sounding over layered structures: Madrid, Interciencia, 150 p., 66 tables. 

Pekeris, C. L., 1940, Direct method of interpretation in resistivity prospecting:

Geophysics, v. 5, no. 1. p. 31-42. 

Rijkswaterstaat, 1969, Standard graphs for resistivity prospecting: The Hague, European

Assoc. Explor. Geophysicists.

Schlichter, L. B., 1933, The interpretation of the resistivity prospecting method for horizontal structures: Physics, v. 4, p. 307-322. 

Shkabarnia, N. G., and Gritsenko, V. G., 1971, Interpretatsiia kriv'ikh elektricheskogo

zondirovaniia s primeneniem EVM: Moscow, Izdatel'ctvo "Nedra," 113 p. 

Stefanesco, S. S., Schlumberger, Conrad, and Schlumberger, Marcel, 1930, Sur la distribution electrique potentielle autor d'une prix de terre ponctuelle dans un terrain a

couche horizontales, homogenes et isotrope: Jour. Physique et Radium, v. 11, no. 1, p.

132-140. 

Strakhov, V. N., 1966a, Numerical methods for the solution of geophysical problems

reducing to integral equations with multiplicative kernels [in Russian]: Akad. Nauk

SSSR Izv. Ser. Geofiz., 1966, no. 9, p. 38-52; translated by F. Goodspeed, in Physics of

the Solid Earth, 1966, p. 570-576. 

___1966b, The analytic determination of the parameters of a horizontally layered

medium from data obtained by vertical electrical prospecting [in Russian]: Akad.

Nauk. SSSR Izv. Ser. Geofiz., 1966, no. 4, p. 52-63; translated by F. Goodspeed, in

Physics of the Solid Earth, 1966 p. 237-242.

.1968, The solution of the inverse problem for vertical electrical sounding [in Russian]: Akad. Nauk. SSSR Izv. Ser. Geofiz., 1968, no. 4, p. 77-84; translated by F.

Goodspeed, in Physics of the Solid Earth, 1968, p. 247-250. 

Strakhov, V. N., and Karelina, G. N., 1969, [On the interpretation of vertical electrical

sounding data with an electronic calculating machine]: Prikladnaya Geofizika, v. 56, p.

118-129 [in Russian]. 

Sunde, E. D., 1949, Earth conduction effects in transmission systems: New York, Van

Nostrand, 370 p. 

Vanyan, L. L., Morozova, G. M., and Lozhenits'na, L., 1962, [On the calculation of

theoretical electrical sounding curves]: Prikladnaya Geofizika, v. 34, p. 135-144 [in

Russian]. 

Vozoff, Keeva, 1958, Numerical resistivity analysis Horizontal layers: Geophysics, v. 23,

no. 3, p. 536-556. 

Watson, G. N., 1962, A treatise on the theory of Bessel functions [2d ed.]: Cambridge,

Cambridge Univ. Press, 804 p. 

Zohdy, A. A. R., 1965, The auxiliary point method of electrical sounding interpretation,

and its relationship to the Dar Zarrouk parameteis: Geophysics, v. 30, no. 4, p. 644-660. 

_____1968, The effect of current leakage and electrode spacing errors on resistivity
INTERPRETATION OF SCHLUMBERGER SOUNDING CURVES E39

measurements, in Geological Survey research 1968: U.S. Geol. Survey Prof. Paper 

600-D, p. D258-D264.

_1974a, Use of Dar Zarrouk curves in the interpretation of vertical electrical

sounding data: U.S. Geol. Survey Bull. 1313-D, 41 p.

_1974b, A computer program for the calculation of Schlumberger sounding curves

by convolution: available only from U.S. Dept. Commerce Natl. Tech. Inf. Service, 

Springfield, Va. 22161, as U.S. Geol. Survey Rept. USGS-GD-74-010, PB-232 056. 

_1974c, A computer program for the automatic interpretation of Schlumberger

sounding curves over horizontally stratified media: available only from U.S. Dept.

Commerce Natl. Tech. Inf. Service, Springfield, Va. 22161, as U.S Geol. Survey Rept.

USGS-GD-74-017, PB-232 703. 

Zohdy, A. A. R., and Jackson, D. B., 1973, Recognition of natural brine by electrical

soundings near the Salt Fork of the Brazos River, Kent and Stonewall Counties, Texas:

U.S. Geol. Survey Prof. Paper 809-A, 14 p. 

Zohdy, A. A. R., and Stanley, W. D., 1973, Preliminary interpretation of electrical

sounding curves obtained across the Snake River Plain from Blackfoot to Arco, Idaho:

U.S. Geol. Survey open-file report, 5 p.

trV.S. GOVERNMENT PRINTING OFFICE: 1983-381-614/586