
Research Article


A Damped LeastSquares Inversion Program for the Interpretation of Schlumberger Sounding Curves 

Yunus Levent Ekinci
and
Alper Demirci



ABSTRACT

An iterative inversion program for 1D interpretation
of Schlumberger sounding curves is presented in this research. The program,
written in MATLAB v.5.3 (R11), optimizes the fit between the observed
and the calculated apparent resistivity data. The procedure is initialized
by userselected initial model parameters and employs the damped leastsquares
solution with the Singular Value Decomposition (SVD). The applicability
of the present inversion program was tested on three synthetic models
which symbolize ideal resistivity distributions of geological relevance.
Noise was added to two models as a further concession to reality. Moreover,
the program was used to interpret real data from Turkey and the results
obtained were compared with published results. In all cases, applications
indicated that the model parameters obtained are in good agreement with
the real ones and the published results. The program is suitable for shallow
small scale investigations if the lateral changes in the study area are
expected to be smooth.







INTRODUCTION
The Direct Current (DC) resistivity sounding method is one of the oldest
geophysical techniques which is intensively used for the evaluation of
the deep and shallow structure of the subsurface. By introducing the electrical
current directly into the ground through a pair of electrodes (current
electrodes), the difference of the resulting voltage can be measured between
the other pair of electrodes (potential electrodes). The apparent resistivity
of the subsurface which is the reciprocal of conductivity can be calculated
by this way. The method aims to determine the resistivity variation with
depth. The depth of the penetration depends on the distance between the
current electrodes. Increasing the depth of the penetration can be carried
out by enlarging the distance between the current electrodes from a small
distance in the beginning to larger distances at the end of the array.
The investigations performed by using DC resistivity sounding method
can yield satisfactory results in the areas where the ground is horizontally
layered with a little lateral variation. The method has been applied to
a wide variety of geophysical exploration problems so far. The most successful
and beneficial surveys have been carried out in hydrogeological, geothermal,
environmental and engineering applications (e.g., Batayneh et al.,
2001; Mota et al., 2004; Agnesi et al., 2005; Özurlan
and Şahin, 2006; Asfahani and Radwan, 2007; Park et al., 2007).
The apparent resistivity data do not represent the true resistivity distribution
of the subsurface. Thus, to determine the true resistivity distribution
of the subsurface, inversion of the observed apparent resistivity data
must be carried out. DC resistivity inverse problem was investigated for
the first time in the early 1930s. From that time until the late 1980s,
the methodology of field survey and character of data from the measurements
remained unchanged. From the late 1980s and early 1990s until today, there
has been considerable improvement in data collection and interpretation
(Muiuane and Pedersen, 1999). In the last decade, depending on high technological
development of computer controlled multielectrode survey systems and
the inversion software packages, much of the investigations carried out
on DC resistivity method have mainly focused on 2D and 3D data acquisition
and interpretation techniques (e.g., Nassir et al., 2000; Dahlin,
2001; Bauer et al., 2006; Drahor et al., 2007; Ekinci et
al., 2008). However, 1D techniques can be used to interpret the data
and may yield good approximation of the earth`s subsurface if the investigation
is carried out in areas with gently sloping or flat lying ground and where
the resistivity and thickness variations between the adjacent layers are
smooth.
In the frame of this study we present a MATLABbased inversion program
that uses damped leastsquares solution for the interpretation of Schlumberger
sounding curves. The reliability of the present program was tested on
both noisy and noisefree synthetic data. Additionally, a field survey
was also conducted in GökçeadaTurkey to scrutinize the efficiency
of the program on real data. We compared the results with our earlier
study which indicates almost horizontally layered resistivity variations
obtained by the interpretation of 2D resistivity inverse model sections
(Ekinci et al., 2008).
MATERIALS AND METHODS
Schlumberger electrode configuration is very simple to use in the survey
area and hence is the most popular DC resistivity sounding array. The
relationship between the apparent resistivity data and the layer parameters
is expressed by an integral equation. The very wellknown expression of
the computation of a forward response for Schlumberger electrode configuration
(ρ_{a}) over an earth model consisting of homogeneous and
isotropic layers is related to the kernel function through the Hankel
integral and can be formulated as follows (Koefoed, 1970):
where, s is half the current electrode spacing in Schlumberger electrode
configuration, J_{1} denotes the firstorder Bessel function of
the first kind and λ denotes the integral variable. The resistivity
transform function, T(λ), is given by the recurrence relationship
as follows:
where, n denotes the number of layers, ρ_{i} and h_{i}
are the resistivity and thickness of the ith layer, respectively. Equation
2 is used to compute the forward model response for DC resistivity
sounding data.
Nonlinear leastsquares inversion scheme iteratively updates the model
parameters in each step with the use of a correction vector which is the
solution of a set of normal equations. Inversion of geoelectrical data
is an illposed problem. Therefore, small changes in the data may lead
large changes in the model. This and the ensuing suboptimality restrict
the initial model to being in the near vicinity of the true model. The
problem may be reduced by introducing damping into the system of equations
(Roy, 1999). This parameter is added to the principal diagonal of A^{T}A
which helps to increase the DC level of eigenvalues such that none of
its eigenvalues can become zero (Raju, 2003). The damped leastsquares
solution is given by the Eq. 3 as follows:
where, Δp is the parameter correction vector, Δd is the data
difference vector, A is the Jacobian matrix. I is the identity matrix
and the term ε is called damping factor (Levenberg, 1944; Marquardt,
1963).
Singular Value Decomposition (SVD) is wellknown technique used in many
areas of applied sciences including the earth sciences. It can be easily
applied to small scale geophysical problems. It is mathematically robust
and numerically stable and also provides other vital information on the
state of the model and data thus enabling model resolution and covariance
studies (Meju, 1994). However it must be noted that the use of SVD is
not logical for large scale problems. The large scale problems can be
solved either explicitly or implicitly using iterative methods like conjugate
gradients. In this study, because the examples used are small scale problems,
Eq. 3 was solved using SVD in the inversion scheme.
An nxn or nxp matrix A can be factored into a product of three other
matrices as follows:
where, for n data and p parameters, U (nxp) and V (pxp) are, respectively
the data space and parameter space eigenvectors and S is a pxp diagonal
matrix containing at most r nonzero eigenvalues of A, with r≤p. These
diagonal entries in S (λ_{1}, λ_{2}, ..., λ_{p})
are called singular values of A. This factorization is known as SVD of
A (Meju, 1994). If SVD is used in damped leastsquares solution, Eq.
5 is obtained:
By adding the damping factor to the diagonal elements we obtain:
The inverse of the Eq. 6 is given by as follows:
Substituting the Eq. 7 in 5 we obtain:
and the parameter correction vector can be expressed as:
Equation 9 provides damped leastsquares solution with
SVD. For automated inversion, the common practice is to set the damping
factor first to a large positive value thus taking the advantage of the
good initial convergence properties of the steepest descent method and
thereafter the damping factor is multiplied by a factor less that unity
after each iteration so that the leastsquares method predominates near
the solution (Meju, 1994). A way of determining the damping factor has
been developed by Arnason and Hersir (1988) and given by as follows:
where, L is the test number performed for the damping factor at any iteration,
λ is the parameter eigenvalue and the term of Δx is given by:
where, x_{r1} is the misfit value obtained previous iteration
and x_{r} is the misfit value of the current iteration. In this
study, Eq. 10 and 11 were used to
set the damping factor in each iteration.
PROGRAM DESCRIPTION
The inversion program presented in this paper is composed of three parts.
One is the main program and the others are the functions.
• 
The main program VES1dinv reads the observed apparent
resistivity data and AB/2 spacings from the data file, uses userselected
initial model parameters (number of layers, resistivity and the thicknesses
of these layers) and inverts the apparent resistivity data by using
damped leastsquares solution with SVD. The program also uses the
basic plotting commands for the illustration of the inversion results. 
• 
The function Jacobian calculates the partial derivatives of parameters
numerically by using a userselected difference interval value and
forms the Jacobian matrix. 
• 
The function VES1dmod calculates the apparent resistivity data from
the given values (resistivity, thickness and AB/2) for Schlumberger
electrode configuration by using the filter coefficients (Nyman and
Landisman, 1977). 
The full codes of the main program and the functions are given in the
Appendix AC, respectively.
SYNTHETIC EXAMPLES
Synthetic models with ideal resistivity distributions of geological relevance
were used in the applications. Three different examples consisting of
three or fourlayered earth models were preferred. In general, as a similar
explanation mentioned by Muiuane and Pedersen (1999), the subsurface models
where resistivity increases or decreases with depth are not difficult
to resolve. On the other hand, if a resistor layer is embedded between
two conductor layers or a conductor layer is embedded between two resistor
layers convergence is more difficult and slow in comparison to first type.
We used these types of models for the numerical examples. The synthetic
data were generated at a predetermined set of 20 AB/2 values for all models
by using the forward modeling function of the program. The 5% noise level
is reasonable for sites with high resistivity contrasts and higher electrode
contact resistances (Dahlin, 1993). Thus, as a further concession to reality,
Model 2 and 3 were contaminated with ±5% random noise depending
on the magnitude of each data point.
First example demonstrates a noisefree threelayered earth model with
decreasing resistivity with depth. Second example deals with a threelayered
earth model consisting of a resistor layer between two conductors. A fourlayered,
decreasing resistivity with depth model was used in the last example.
We used homogeneous half spaces for the initial models. The inversion
procedure was controlled by two misfit functions in the iterations as
given in Appendix A. When the predetermined error limit was reached, the
inversion procedure was stopped automatically. The convergences were calculated
as an rms error function between the observed and calculated apparent
resistivity values at the end of the process as follows:
where, obs is observed, cal is calculated apparent resistivity values
and ND is the number of the apparent resistivity data. In all cases, although
the initial guesses were far removed from the true ones, the model parameters
were resolved and the program yielded satisfactory results. The inversion
results are illustrated as bestfit curves and resistivitydepth distributions
for Model 1, 2 and 3 in Fig. 13, respectively.

Fig. 1: 
Model 1, (a) observed and calculated apparent resistivity
curves and (b) resistivitydepth distribution of the initial and final
models 

Fig. 2: 
Model 2, (a) observed and calculated apparent resistivity
curves and (b) resistivitydepth distribution of the initial and final
models 

Fig. 3: 
Model 3, (a) observed and calculated apparent resistivity
curves and (b) resistivitydepth distribution of the initial and final
models 
Table 1: 
Observed and calculated apparent resistivity data of
the models 

Table 2: 
True and obtained model parameters of the models 

Observed
and calculated apparent resistivity values and obtained model parameters
after inversion procedure are given in Table 12, respectively.
REAL DATA EXAMPLE
To test the validity of the present inversion program on real data, a
Schlumberger sounding curve was interpreted. The survey area, named Kaleköy
is located in the northeastern part of Gökçeada (Fig.
4). The area has been previously studied (Ekinci et al., 2008).
The seawater and the freshwater saturated layer were determined by the
interpretation of 2D electrical resistivity tomography survey. The resistivity
model section shows almost horizontally layered resistivity variations
(Ekinci et al., 2008). Therefore, we expected that Schlumberger
sounding should provide similar results.
The geological settings of Gökçeada are represented by a
thick sedimentary sequence, volcanic rock and alluvium. Regarding of the
survey site, western and eastern parts of it in which Mezardere Formation
outcrops are hilly and the formation consists of sandstoneshale sequence.
The plain is also represented by alluvium (Fig. 4).
Hill slope allows rapid surface runoff wherefore infiltration occurs easily
into alluvium. The existence of water saturation in this area is suggested
due to both hydrogeological conditions (shallow water level, rapid surface
runoff) and the high Pwave velocity in the alluvium (Ekinci et al.,
2008). The site is also appropriate to investigate the seawater existence.
Thus, we applied a Schlumberger sounding in order to test the efficiency
of the developed program on the determination of seawater existence and
its depth and thickness. The measurements were performed with Iris SyscalR1+
resistivity meter. The maximum halfcurrent electrode separation was 120
m. In order to improve the quality of the data, the resistivity meter
performed a stack of a minimum of three for each data point.

Fig. 4: 
The location map of site and the Schlumberger sounding
(not to scale) 

Fig. 5: 
(a) observed and calculated apparent resistivity curves
and (b) resistivitydepth distribution of the initial and final models
for the field example 
Table 3: 
The obtained model parameters of the real data 

When the
standard deviation between the measurements reached 3%, a maximum of five
measurements were performed for each data point during the investigation.
The standard deviations were generally below 1%.
We used a fourlayered earth model for the initial guess. The inversion
converged after the 65 th iteration with an rms error of 0.77. The fit
between the observed and calculated apparent resistivity data are shown
in Fig. 5 and obtained model parameters are given in
Table 3. The sounding curve represents decreasing resistivity
with depth but ends in a resistive basal. According to the geological
evaluation of the inversion results, the top layer represents the consolidated
alluvium of ~2.5 m overlying a water saturated zone, which extends to
a depth of ~13.5 m. This water saturated zone exists above the very conductive
layer (< 2.5 Ωm). The source of the freshwater is thought to be
the groundwater recharge process, which occurs directly from rainfall
and from rapid surface runoff wherefore infiltration occurs easily into
alluvium due to the hill slopes. The very conductive layer which extends
to a depth of ~23 m corresponds to seawater. Additionally, the increase
of resistivity values below that layer may trace the existence of Mezardere
formation. This may be interpreted to be decreasing porosity in response
to increasing compaction and this compaction may indicate Mezardere formation
represented by sandstoneshale sequence. These results are in close agreement
with those obtained by Ekinci et al. (2008).
CONCLUSIONS
A MATLABbased 1D inversion program for the interpretation of the Schlumberger
sounding curves was presented with synthetic and real data examples in
this paper. The presented inversion scheme is computationally simple,
fast and easy to carry out. The program uses damped leastsquares solution
with SVD and runs iteratively by updating the model parameters.
The program was tested on small scale synthetically produced noisefree
and noisy apparent resistivity data sets which have different curve types.
Although, the initial guesses were far removed from the true ones, the
present inversion program produced satisfactory results. Additionally,
the results obtained by the inversion of the real data are in close agreement
with the results of an earlier study. The existence of seawater and its
depth and thickness were easily determined.
The program is suitable to use by the interpreter for small scale problems
where the field conditions do not allow gathering enough data for multidimensional
modeling scheme. It has been concluded that the presented program can
provide reliable results where the resistivity and thickness variations
between the adjacent layers are smooth. This inversion program is developed
for Schlumberger electrode configuration but due to the fact that it is
an open source code it can be easily modified for the other electrode
configurations by changing the forward modeling function or extended to
the needs of the other users.
ACKNOWLEDGMENTS
The authors would like to express their gratitude to Dr. E.U. Ulugergerli for
his useful suggestions. Thanks to Dr. T. Bekler for his permission to use the
data set. They also thank to C. Ertekin, scientist, for his valuable information
about the geological and hydrogeological properties of the surveyed area. The
following undergraduate students are thanked for their aid in the acquisition
of the field data: H. Şahlan, M. Özben, S. Kutluk, S. Güneş
and O. Zop. The authors are especially grateful to Dr. M.A. Riahi for his critiques
that have helped us to improve the paper and also thanks to three anonymous
reviewers for their valuable comments. APPENDIX A
% VES1dinv % 
% 1D Inversion of Schlumberger Sounding Data % 
% By Yunus Levent Ekinci and Alper Demirci % 
close all; 
clear all; 
clc; 
format long; 
load ab2.dat; 
x = ab2(:,1); 
roa = ab2(:,2); 
r = [ρ_{1} ρ_{2} ρ_{n}]; 
t = [t_{1} t_{n1}]; 
m = [r t]; 
rinitial = r; 
tinitial = t; 
lr = length(r); 
lt = length(t); 
kr = 10e10; 
iteration = 1; 
maxiteration = 500; 
dfit = 1; 
while iteration<maxiteration 
r = m(1:lr); 
t = m(1+lr:lr+lt); 
for i = 1:length(x) 
s = ab2(i); 
[g] = VES1dmod (r,t,s); 
roa1(i,:) = g; 
end 
e1 = [log(roa)log(roa1)]; 
dd = e1; 
misfit1 = e1`*e1; 
if misfit1<kr 
loglog(x,roa,`k.`,x,roa1,`k`); 
axis([1 1000 1 1000]) 
xlabel(`AB/2 (m)`); 
ylabel(`Apparent Resistivity (Ohmm)`); 
break 
end 
[A] = jacobian(ab2,x,r,t,lr,lt,roa,roa1); 
[U S V] = svd(A,0); 
ss = length(S); 
say = 1; 
k = 0; 
while say<ss 
diagS = diag(S); 
beta = S(say)*(dfit^(1/say)); 
if beta<10e5 
beta = 0.001*say; 
end 
for i4 = 1:ss 
SS(i4,i4) = S(i4,i4)/(S(i4,i4)^2+beta); 
end 
dmg = V*SS*U`*dd; 
mg = exp(log(m)+dmg`); 
r = mg(1:lr); 
t = mg(1+lr:lr+lt); 
for i5 = 1:length(x) 
s = ab2(i5); 
[g] = VES1dmod (r,t,s); 
roa4(i5,:) = g; 
end 
e2 = [log(roa)log(roa4)]; 
misfit2 = e2`*e2; 
if misfit2>misfit1 
(`Beta control`) 
say = say+1; 
k = k+1; 
if k = = ss1 
iteration = maxiteration; 
say = ss+1; 
end 
else 
say = ss+1; 
m = mg; 
dfit = (misfit1misfit2)/misfit1; 
iteration = iteration+1; 
a = iteration; 
if dfit<kr 
iteration = maxiteration; 
say = say+1; 
end 
end 
end 
subplot(1,2,1), 
loglog(x,roa,`k.`,x,roa4,`k`); 
axis([1 1000 1 1000]) 
xlabel(`AB/2 (m)`); 
ylabel(`Apparent Resistivity (Ohmm)`); 
pause(0.001) 
end 
observed = roa; 
calculated = roa4; 
h = legend(`obs`,`cal`,1); 
format bank; 
m 
rr = [0,r]; 
tt = [0,cumsum(t),max(t)*10]; 
subplot(1,2,2), 
stairs(rr,tt,`k`); 
rrr = [0,rinitial]; 
ttt = [0,cumsum(tinitial),max(t)*10]; 
hold on; 
subplot(1,2,2), 
stairs(rrr,ttt,`k`); 
set(gca,`Ydir`,`reverse`); 
set(gca,`Xscale`,`log`); 
ylim([0 50]); 
xlim([0 1000]); 
xlabel(`Resistivity (Ohmm)`); 
ylabel(`Depth (m)`); 
rms = norm(roa4roa) /sqrt(length(roa)); 
APPENDIX B
% Function Jacobian % 
function[A] = jacobian(ab2,x,r,t,lr,lt,roa,roa1) 
par = 0.1; 
r2 = r; 
for i2 = 1:lr 
r2(i2) = (r(i2)*par)+r(i2); 
for ii = 1:length(x) 
s = ab2(ii); 
[g] = VES1dmod (r2,t,s); 
roa2(ii,:) = g; 
end 
A1(:,i2) = [(roa2roa1)/(r(i2)*par)]*r(i2)./roa; 
r2 = r; 
end 
t2 = t; 
for i3 = 1:lt 
t2(i3) = (t(i3)*par)+t(i3); 
for ii = 1:length(x) 
s = ab2(ii); 
[g] = VES1dmod (r,t2,s); 
roa3(ii,:) = g; 
end 
A2(:,i3) = [(roa3roa1)/(t(i3)*par)]*t(i3)./roa; 
t2 = t; 
end 
A = [A1 A2]; 
return 
APPENDIX C
% Function VES1dmod % 
function[g] = VES1dmod(r,t,s) 
q = 13; 
f = 10; 
m = 4.438; 
x = 0; 
e = exp(0.5*log(10)/m); 
h = 2*q2; 
u = s*exp(f*log(10)/mx); 
l = length(r); 
n = 1; 
for i = 1:n+h 
w = l; 
v = r(l); 
while w>1 
w = w1; 
aa = tanh(t(w)/u); 
v = (v+r(w)*aa)/(1+v*aa/r(w)); 
end 
a(i) = v; 
u = u*e; 
end 
i = 1; 
g = 105*a(i)262*a(i+2)+416*a(i+4)746*a(i+6)+1605*a(i+8); 
g = g4390*a(i+10)+13396*a(i+12)27841*a(i+14); 
g = g+16448*a(i+16)+8183*a(i+18)+2525*a(i+20); 
g = (g+336*a(i+22)+225*a(i+24))/10000; 
return 

REFERENCES 
1: Agnesi, V., M. Camarda, C. Conoscenti, C. Di Maggio, I.S. Diliberto, P. Madonia and E. Rotigliano, 2005. A multidisciplinary approach to the evaluation of the mechanism that triggered the Cerda landslide (Sicily, Italy). Geomorphology, 65: 101116. CrossRef  Direct Link 
2: Asfahani, J. and Y. Radwan, 2007. Tectonic evolution and hydrogeological characteristic of the Khanaser Valley, Northern Syria, derived from the interpretation of vertical electrical sounding. P. Applied Geophys., 164: 22912311. CrossRef  Direct Link 
3: Batayneh, A.T., A.S. Al Zoubi and A.A. Abueladas, 2001. Geophysical investigations for the location of a proposed dam in Al Bishriyya (Al Aritayn) area, Northeast Badia of Jordan. Environ. Geol., 40: 918922. CrossRef  Direct Link 
4: Bauer, P., R. Supper, S. Zimmermann and W. Kinzelbach, 2006. Geoelectrical imaging of groundwater salinization in the Okavango Delta, Botswana. J. Applied Geophys., 60: 126141. CrossRef  Direct Link 
5: Dahlin, T., 1993. On the automation of 2D resistivity surveying for engineering and environmental applications. Ph.D Thesis, Lund Institut of Technology.
6: Dahlin, T., 2001. The development of DC resistivity imaging techniques. Comput. Geosci., 27: 10191029. CrossRef  Direct Link 
7: Drahor, M.G., G. Gokturkler, M.A. Berge, T.O. Kutulmq and N. Tuna, 2007. 3D resistivity imaging from an archaeological site in SouthWestern Anatolia: A case study. Near Surf. Geophys., 5: 195201. Direct Link 
8: Ekinci, Y.L., A. Demirci and C. Ertekin, 2008. Delineation of the seawaterfreshwater interface from the coastal alluvium of KalekoyGokceada, NW Turkey. J. Applied Sci., 8: 19771981. CrossRef  Direct Link 
9: Koefoed, O., 1970. A fast method for determining the layer distribution from the raised kernel function in geoelectrical soundings. Geophys. Prospect., 18: 564570. CrossRef  Direct Link 
10: Levenberg, K., 1944. A method for the solution of certain nonlinear problems in leastsquares. Q. Applied Math., 2: 164168. Direct Link 
11: Marquardt, D.W., 1963. An algorithm for leastsquares estimation of non linear parameters. J. Soc. Ind. Applied Math., 11: 431441. Direct Link 
12: Nyman, D.C. and M. Landisman, 1977. VES dipoledipole filter coefficients. Geophysics, 42: 10371044. CrossRef  Direct Link 
13: Meju, M.A., 1994. Geophysical Data Analysis: Understanding Inverse Problem Theory and Practice. Society of Exploration Geophysics Course Notes Series, No. 6, 1st Edn., SEG Publishers, Tulsa, Oklahoma, ISBN: 1560800275, pp: 296.
14: Mota, R., F.A.M. Santos, A. Mateus, F.Q. Marques, M.A. Goncalves, J. Figueiras and H. Amaral, 2004. Granite fracturing and incipient pollution beneath a recent landfill facility as detected by geoelectrical surveys. J. Applied Geophys., 57: 1122. CrossRef  Direct Link 
15: Muiuane, E.A. and L.B. Pedersen, 1999. Automatic 1D interpretation of DC resistivity sounding data. J. Applied Geophys., 42: 3545. CrossRef  Direct Link 
16: Nassir, S.S.A., M.H. Loke, C.Y. Lee and M.N.M. Nawawi, 2000. Saltwater intrusion mapping by geoelectrical imaging surveys. Geophys. Prospect., 48: 647661. CrossRef  Direct Link 
17: Ozurlan, G. and M.H. Sahin, 2006. Integrated geophysical investigations in the Hisar geothermal field, Demirci, Western Turkey. Geothermics, 35: 110122. CrossRef  Direct Link 
18: Park, Y.H., S.J. Doh and S.T. Yun, 2007. Geoelectric resistivity sounding of riverside alluvial aquifer in an agricultural area at Buyeo, Geum River watershed, Korea: An application to groundwater contamination study. Environ. Geol., 53: 849859. CrossRef  Direct Link 
19: Raju, D.C.V., 2003. LIMAT: A computer program for leastsquares inversion of magnetic anomalies over long tabular bodies. Comput. Geosci., 29: 9198. CrossRef  Direct Link 
20: Roy, I.G., 1999. An efficient nonlinear leastsquares 1D inversion scheme for resistivity and IP sounding data. Geophys. Prospect., 47: 527550. CrossRef  Direct Link 
21: Arnason, K. and G.P. Hersir, 1988. One dimensional inversion of Schlumberger resistivity soundings: Computer Program, Description and Userâ€™s Guide: The United Nations University, Geothermal Training, Report 8, pp: 59.



