2-D and 3-D Electrical Imaging Surveys
2-D and 3-D Electrical Imaging Surveys
2-D and 3-D Electrical Imaging Surveys
By
Dr. M.H.Loke
Copyright (1996-2018)
The author, M.H.Loke, retains the copyright to this set of notes. Users may print a copy
of the notes, but may not alter the contents in any way. The copyright notices must be
retained. For public distribution, prior approval by the author is required.
It is hoped that the information provided will prove useful for those carrying out 2-D
and 3-D field surveys, but the author will not assume responsibility for any damage or
loss caused by any errors in the information provided. If you find any errors, please
inform me by email and I will make every effort to correct it in the next edition.
You can download the programs mentioned in the text (RES2DMOD, RES2DINV,
RES3DMOD, RES3DINV) from the following Web site
www.geotomosoft.com
M.H.Loke
Table of Contents
1 Introduction to resistivity surveys ....................................................................................... 1
1.1 Introduction and basic resistivity theory ...................................................................... 1
1.2 Electrical properties of earth materials ......................................................................... 5
1.3 1-D resistivity surveys and inversions – applications, limitations and pitfalls ............ 6
1.4 Basic Inverse Theory .................................................................................................. 11
1.5 2-D model discretization methods .............................................................................. 15
2 2-D electrical surveys – Data acquisition, presentation and arrays ................................... 18
2.1 Introduction ................................................................................................................ 18
2.2 Field survey method - instrumentation and measurement procedure ......................... 18
2.3 Available field instruments ......................................................................................... 21
2.4 Pseudosection data plotting method ........................................................................... 24
2.5 A comparison of the different electrode arrays .......................................................... 26
2.5.1 The Frechet derivative for a homogeneous half-space ........................................ 26
2.5.2 A 1-D view of the sensitivity function - depth of investigation .......................... 27
2.5.3 A 2-D view of the sensitivity function ................................................................. 30
2.5.4 Wenner array........................................................................................................ 31
2.5.5 Dipole-dipole array .............................................................................................. 31
2.5.6 Wenner-Schlumberger array ................................................................................ 35
2.5.7 Pole-pole array ..................................................................................................... 37
2.5.8 Pole-dipole array .................................................................................................. 38
2.5.9 Multiple gradient array ........................................................................................ 41
2.5.10 High-resolution electrical surveys with overlapping data levels ...................... 44
2.5.11 Summary of array types .................................................................................... 45
2.5.12 Use of the sensitivity values for multi-channel measurements or streamers .... 46
3 A 2-D forward modeling program .................................................................................... 48
3.1 Finite-difference and finite-element methods............................................................. 48
3.2 Using the forward modeling program RES2DMOD .................................................. 49
3.3 Forward modeling exercises ....................................................................................... 50
4 A 2-D inversion program .................................................................................................. 52
4.1 Introduction ................................................................................................................ 52
4.2 Pre-inversion and post-inversion methods to remove bad data points ....................... 52
4.3 Selecting the proper inversion settings ....................................................................... 55
4.4 Using the model sensitivity and uncertainty values ................................................ 61
4.5 Methods to handle topography ................................................................................... 66
4.6 Incorporating information from borehole logs and seismic surveys .......................... 68
4.7 Model refinement ....................................................................................................... 72
4.8 Fast inversion of long 2-D survey lines ...................................................................... 76
4.8.1 Preprocessing steps .............................................................................................. 76
4.8.2 Data set and computer system used for tests ....................................................... 76
4.8.3 Finite-element mesh size...................................................................................... 77
4.8.4 Limit calculation of Jacobian matrix values ........................................................ 77
4.8.5 Use sparse inversion techniques .......................................................................... 78
4.8.6 Use wider cells ..................................................................................................... 78
4.8.7 Overview of methods to reduce the calculation time........................................... 79
4.9 Model resolution and automatic array optimization methods .................................... 81
4.9.1 Concept of model resolution ................................................................................ 81
4.9.2 A heuristic explanation of model resolution – through a glass darkly ................ 81
4.9.3 Examples of model resolution with standard arrays ............................................ 82
4.9.4 Array optimization methods ................................................................................ 83
4.9.5 DOI versus model resolution? ............................................................................. 86
4.9.6 The streamer design problem using a model resolution approach. ...................... 87
4.10 Pitfalls in 2-D resistivity surveys and inversion...................................................... 89
4.11 The pole-pole inversion paradox ............................................................................. 94
5 I.P. inversion ..................................................................................................................... 96
5.1 Introduction ................................................................................................................ 96
5.2 The IP effect ............................................................................................................... 96
5.3 IP data types................................................................................................................ 98
5.4 IP mathematical models............................................................................................ 100
5.5 I.P. surveys with multi-electrode systems ................................................................ 101
6 Cross-borehole imaging .................................................................................................. 103
6.1 Introduction .............................................................................................................. 103
6.2 Electrode configurations for cross-borehole surveys ............................................... 103
6.2.1 Two electrodes array – the pole-pole ................................................................. 103
6.2.2 Three electrodes array – the pole-bipole ............................................................ 105
6.2.3 Four electrodes array – the bipole-bipole .......................................................... 107
6.3 Single borehole surveys ............................................................................................ 110
6.4 Cross-borehole optimized arrays .............................................................................. 112
6.5 Optimized arrays with subsurface electrodes ........................................................... 112
7 2-D field examples .......................................................................................................... 116
7.1 Introduction .............................................................................................................. 116
7.2 Landslide - Cangkat Jering, Malaysia ...................................................................... 116
7.3 Old Tar Works - U.K. ............................................................................................... 116
7.4 Holes in clay layer - U.S.A. ...................................................................................... 117
7.5 Time-lapse water infiltration survey - U.K. ............................................................. 119
7.6 Pumping test, U.K. ................................................................................................... 121
7.7 Wenner Gamma array survey - Nigeria .................................................................... 122
7.8 Mobile underwater survey - Belgium ....................................................................... 123
7.9 Floating electrodes survey – U.S.A. ......................................................................... 125
7.10 Oil Sands, Canada ................................................................................................. 125
8 3-D electrical imaging surveys ........................................................................................ 127
8.1 Introduction to 3-D surveys ...................................................................................... 127
8.2 Array types for 3-D surveys ..................................................................................... 127
8.2.1 The pole-pole array ............................................................................................ 127
8.2.2 The pole-dipole array ......................................................................................... 130
8.2.3 The dipole-dipole, Wenner and Schlumberger arrays ....................................... 130
8.2.4 Summary of array types ..................................................................................... 132
8.3 3-D roll-along techniques ......................................................................................... 135
8.4 A 3-D forward modeling program ............................................................................ 137
8.5 3-D inversion algorithms and 3-D data sets ............................................................. 139
8.6 A 3-D inversion program .......................................................................................... 140
8.7 Banding effects in 3-D inversion models ................................................................. 141
8.8 The use of long electrodes in 3-D surveys ............................................................... 144
8.9 Data grid formats and model discretizations ............................................................ 146
8.9.1 Types of surveys and model grids ..................................................................... 146
8.9.2 Model grid optimization for large surveys with arbitrary electrodes positions . 147
8.10 3-D array optimization – grids and perimeters...................................................... 151
8.11 Unstable arrays and the geometric factor relative error ........................................ 158
8.12 Not on firm foundations – inversion with shifting electrodes in 2-D and 3-D ..... 160
8.13 Examples of 3-D field surveys .............................................................................. 165
8.13.1 Birmingham field test survey - U.K. .............................................................. 165
8.13.2 Sludge deposit - Sweden................................................................................. 166
8.13.4 Copper Hill - Australia ................................................................................... 169
8.13.5 Athabasca Basin – Canada ............................................................................. 171
8.13.6 Cripple Creek and Victor Gold Mine, Colorado – USA ................................ 173
8.13.7 Burra copper deposit, South Australia : Model reliability determination ...... 175
8.14 Closing remarks on the 3-D method ..................................................................... 178
Acknowledgments.................................................................................................................. 179
References .............................................................................................................................. 180
Appendix A The smoothness constraint and resolving deeper structures in I.P. surveys 189
A.1 A problem with I.P. inversion with a conductive overburden .................................. 189
A.2 2-D synthetic model test ........................................................................................... 189
A.3 3-D field data set test ................................................................................................ 190
A.4 Software implementation.......................................................................................... 192
Appendix B Modeling long electrodes ............................................................................ 193
B.1 Two methods to model a long electrode ................................................................... 193
Appendix C New I.P. systems and negative apparent resistivity values ......................... 194
C.1 The distributed receivers I.P. systems and negative apparent resistivity values. ..... 194
List of Figures
Figure 1. The flow of current from a point current source and the resulting potential
distribution. ......................................................................................................................... 3
Figure 2. The potential distribution caused by a pair of current electrodes. The electrodes are
1 m apart with a current of 1 ampere and a homogeneous half-space with resistivity of 1
m...................................................................................................................................... 3
Figure 3. A conventional array with four electrodes to measure the subsurface resistivity. .... 4
Figure 4. Common arrays used in resistivity surveys and their geometric factors. Note that
the dipole-dipole, pole-dipole and Wenner-Schlumberger arrays have two parameters, the
dipole length “a” and the dipole separation factor “n”. While the “n” factor is commonly
an integer value, non-integer values can also be used......................................................... 4
Figure 5. The resistivity of rocks, soils and minerals. ............................................................... 6
Figure 6. The three different models used in the interpretation of resistivity measurements. ... 7
Figure 7. A typical 1-D model used in the interpretation of resistivity sounding data for the
Wenner array. ...................................................................................................................... 7
Figure 8. A 2-D two-layer model with a low resistivity prism in the upper layer. The
calculated apparent resistivity pseudosections for the (a) Wenner and (b) Schlumberger
arrays. (c) The 2D model. The mid-point for a conventional sounding survey is also
shown. ................................................................................................................................. 9
Figure 9. Apparent resistivity sounding curves for a 2-D model with a lateral inhomogeneity.
(a) The apparent resistivity curve extracted from the 2D pseudosection for the Wenner
array. The sounding curve for a two-layer model without the low resistivity prism is also
shown by the black line curve. (b) The apparent resistivity curves extracted from the 2-D
pseudosection for the Schlumberger array with a spacing of 1.0 meter (black crosses) and
3.0 meters (red crosses) between the potential electrodes. The sounding curve for a two-
layer model without the low resistivity prism is also shown. ........................................... 10
Figure 10. Coupling between neighboring model cells through the roughness filter in a 2-D
model. (a) In the horizontal and vertical diretcions only, and (b) in diagonal diretcions as
well. ................................................................................................................................... 13
Figure 11. (a) Coupling between corresponding model blocks in two time-lapse models using
a cross-model time-lapse smoothness constraint. (b) Example Jacobian matrix structure
for five time series data sets and models. Each grey rectangle represents the Jacobian
matrix associated with a single set of measurements. ....................................................... 15
Figure 12. The different models for the subsurface used in the interpretation of data from 2-D
electrical imaging surveys. (a) A purely cell based model. (b) A purely boundary based
model. (c) The laterally constrained model. (d) A combined cell based and boundary
based model with rectangular cells, and (e) with boundary conforming trapezoidal cells.
........................................................................................................................................... 16
Figure 13. Example of a cell based model with a variable boundary. (a) The test model. (b)
The apparent resistivity pseudosection. (c) The inversion model with a variable sharp
boundary that is marked by a black line. ........................................................................... 17
Figure 14. The arrangement of electrodes for a 2-D electrical survey and the sequence of
measurements used to build up a pseudosection. .............................................................. 19
Figure 15. The use of the roll-along method to extend the area covered by a 2-D survey. ..... 20
Figure 16. Sketch outline of the ABEM Lund Imaging System. Each mark on the cables
indicates an electrode position (Dahlin, 1996). The cables are placed along a single line
(the sideways shift in the figure is only for clarity). This figure also shows the principle
of moving cables when using the roll-along technique. .................................................... 22
Figure 17. The Aarhus Pulled Array System. The system shown has two current (C)
electrodes and six potential electrodes (Christensen and Sørensen 1998, Bernstone and
Dahlin 1999)...................................................................................................................... 23
Figure 18. The Geometrics OhmMapper system using capacitive coupled electrodes.
(Courtesy of Geometrics Inc.). .......................................................................................... 23
Figure 19. Schematic diagram of a possible mobile underwater survey system. The cable has
two fixed current electrodes and a number of potential electrodes so that measurements
can be made at different spacings. The above arrangement uses the Wenner-
Schlumberger type of configuration. Other configurations, such as the gradient array, can
also be used. ...................................................................................................................... 24
Figure 20. The apparent resistivity pseudosections from 2-D imaging surveys with different
arrays over a rectangular prism. ........................................................................................ 25
Figure 21. The parameters for the sensitivity function calculation at a point (x,y,z) within a
half-space. A pole-pole array with the current electrode at the origin and the potential
electrode “a” meters away is shown. ................................................................................ 26
Figure 22. A plot of the 1-D sensitivity function. (a) The sensitivity function for the pole-pole
array. Note that the median depth of investigation (red arrow) is more than twice the
depth of maximum sensitivity (blue arrow). (b) The sensitivity function and median
depth of investigation for the Wenner array...................................................................... 28
Figure 23. 2-D sensitivity sections for the Wenner array. The sensitivity sections for the (a)
alpha, (b) beta and (c) gamma configurations. .................................................................. 32
Figure 24. 2-D sensitivity sections for the dipole-dipole array. The sections with (a) n=1, (b)
n=2, (c) n=4 and (d) n=6. ................................................................................................. 33
Figure 25. Two possible different arrangements for a dipole-dipole array measurement. The
two arrangements have the same array length but different “a” and “n” factors resulting
in very different signal strengths. ...................................................................................... 35
Figure 26. 2-D sensitivity sections for the Wenner-Schlumberger array. The sensitivity
sections with (a) n=1, (b) n=2, (c) n=4 and (d) n=6. ........................................................ 36
Figure 27. A comparison of the (i) electrode arrangement and (ii) pseudosection data pattern
for the Wenner and Wenner-Schlumberger arrays. ........................................................... 37
Figure 28. The pole-pole array 2-D sensitivity section............................................................ 38
Figure 29. The forward and reverse pole-dipole arrays. .......................................................... 39
Figure 30. The pole-dipole array 2-D sensitivity sections. The sensitivity sections with (a)
n=1, (b) n=2, (c) n=4 and (d) n=6. .................................................................................... 40
Figure 31. Sensitivity sections for the gradient array. The current electrodes are fixed at x=0
and 1.0 meters, and the distance between the potential electrodes is 0.1 m. The position
of the P1-P2 potential electrodes at (a) 0.45 and 0.55 m., (b) 0.55 and 0.65 m., (c) 0.65
and 0.75 m, (d) 0.75 and 0.85 m. and (e) 0.80 and 0.90 m. .............................................. 42
Figure 32. (a) Example of multiple gradient array data set and inversion model. (b) Profile
plot using exact pseudodepths. (c) Profile plot using approximate pseudodepths. ........... 43
Figure 33. The apparent resistivity pseudosection for the dipole–dipole array using
overlapping data levels over a rectangular prism. Values of 1 to 3 meters are used for the
dipole length ‘a’, and the dipole separation factor ‘n’ varies from 1 to 5. Compare this
with Figure 20c for the same model but with ‘a’ fixed at 1 meter, and “n” varying from 1
to 10. In practice, a ‘n’ value greater than 8 would result in very noisy apparent resistivity
values. ................................................................................................................................ 45
Figure 34. Cumulative sensitivity sections for different measurement configurations using (a)
a dipole-dipole sequence, (b) a moving gradient array and (c) an expanding Wenner-
Schlumberger array. .......................................................................................................... 47
Figure 35. The output from the RES2DMOD software for the SINGLE_BLOCK.MOD 2-D
model file. The individual cells in the model are shown in the lower figure, while the
upper figure shows the pseudosection for the Wenner Beta (dipole-dipole with n=1)
array. .................................................................................................................................. 49
Figure 36. An example of a field data set with a few bad data points. The most obvious bad
datum points are located below the 300 meters and 470 meters marks. The apparent
resistivity data in (a) pseudosection form and in (b) profile form. ................................... 53
Figure 37. Selecting the menu option to remove bad data points manually. ........................... 54
Figure 38. Error distribution bar chart from a trial inversion of the Grundfor Line 1 data set
with five bad data points. .................................................................................................. 55
Figure 39. Different options to modify the inversion process. ................................................ 56
Figure 40. Example of inversion results using the l 2-norm smooth inversion and l1-norm
blocky inversion model constrains. (a) Apparent resistivity pseudosection (Wenner array)
for a synthetic test model with a faulted block (100 m) in the bottom-left side and a
small rectangular block (2 m) on the right side with a surrounding medium of 10 m.
The inversion models produced by (b) the conventional least-squares smoothness-
constrained or l2-norm inversion method and (c) the robust or l1-norm inversion method.
........................................................................................................................................... 57
Figure 41. Different methods to subdivide the subsurface into rectangular prisms in a 2-D
model. Models obtained with (a) the default algorithm, (b) by allowing the number of
model cells to exceed the number of data points, (c) a model which extends to the edges
of the survey line and (d) using the sensitivity values for a homogeneous earth model. .. 58
Figure 42. The options to change the thickness of the model layers. ...................................... 59
Figure 43. The options under the ‘Change Settings’ menu selection. ..................................... 60
Figure 44. The dialog box to limit the model resistivity values. ............................................. 61
Figure 45. Landfill survey example (Wenner array). (a) Apparent resistivity pseudosection ,
(b) model section, (c) model sensitivity section (d) model uncertainty section, (e)
minimum and maximum resistivity sections. .................................................................... 62
Figure 46. Landfill survey depth of investigation determination. (a) Model section with
extended depths and (b) the normalized DOI index section. ............................................ 65
Figure 47. Beach survey example in Denmark. The figure shows the model section with
extended depths and the normalized DOI index section. .................................................. 65
Figure 48. Different methods to incorporate topography into a 2-D inversion model. (a)
Schematic diagram of a typical 2-D inversion model with no topography. A finite-
element mesh with four nodes in the horizontal direction between adjacent electrodes is
normally used. The near surface layers are also subdivided vertically by several mesh
lines. Models with a distorted grid to match the actual topography where (b) the
subsurface nodes are shifted vertically by the same amount as the surface nodes, (c) the
shift in the subsurface nodes are gradually reduced with depth or (d) rapidly reduced with
depth, and (e) the model obtained with the inverse Schwartz-Christoffel transformation
method. .............................................................................................................................. 67
Figure 49. Fixing the resistivity of rectangular and triangular regions of the inversion model.
........................................................................................................................................... 70
Figure 50. The inversion model cells with fixed regions. The fixed regions are drawn in
purple. Note that the triangular region extends beyond the survey line. ........................... 71
Figure 51. Example of an inversion model with specified sharp boundaries. (a) The
boundaries in the Clifton survey (Scott et al., 2000) data set is shown by the blue lines.
(b) The measured apparent resistivity pseudosection and the inversion model section. ... 71
Figure 52. The effect of cell size on the model misfit for near surface inhomogeneities. (a)
Model with a cell width of one unit electrode spacing. (b) A finer model with a cell width
of half the unit electrode spacing. The near surface inhomogeneities are represented by
coloured ovals. .................................................................................................................. 72
Figure 53. Synthetic model (c) used to generate test apparent resistivity data for the pole-
dipole (a) and Wenner (b) arrays....................................................................................... 73
Figure 54. The effect of cell size on the pole-dipole array inversion model. Pole-dipole array.
(a) The apparent resistivity pseudosection. The inversion models obtained using cells
widths of (b) one, (b) one-half and (c) one-quarter the unit electrode spacing. ................ 74
Figure 55. Example of the use of narrower model cells with the Wenner-Schlumberger array.
(a) The apparent resistivity pseudosection for the PIPESCHL.DAT data set. The
inversion models using (b) cells with a width of 1.0 meter that is the same as the actual
unit electrode, and (c) using narrower cells with a width of 0.5 meter. ............................ 75
Figure 56. Example of reduction of near-surface 'ripples' in inversion model. (a) Apparent
resistivity pseudosection for the BLUERIDGE.DAT data set. (b) Normal inversion using
the robust inversion norm and model refinement. (c) Inversion using higher damping
factor for first layer to reduce the 'ripple' effect in the top layer. ...................................... 75
Figure 57. Inversion models for a long survey line using different settings to reduce the
calculation time. (a) The measured apparent resistivity pseudosection for the Redas
underwater survey for the first 2600 meters. Inversion models using (b) standard
inversion settings with 4 nodes between adjacent electrode positions, (c) standard
inversion settings with 2 nodes between adjacent electrode positions, (d) with calculation
of Jacobian matrix values for selected model cells, (e) with sparse inversion technique,
(f) with 2 meter model cell width and (g) with 3 meter model cell width. ....................... 80
Figure 58. Part of finite-element mesh used to model a survey with submerged electrodes.
The resistivity of mesh cells in the water layer are fixed at the known water resistivity. 80
Figure 59. Part of a finite-element mesh for a long survey line. The electrodes used in the
same array are 3 meters apart, while the data unit electrode spacing is 1 meter. Using 2
nodes between adjacent electrodes in the mesh will actually result in 6 nodes between
electrodes in the same array. ............................................................................................. 81
Figure 60. Model resolution sections for the (a) Wenner, (b) Wenner-Schlumberger and (c)
simple dipole-dipole array (d) dipole-dipole array with overlapping data levels. ............ 83
Figure 61. (a) Model resolution section for comprehensive data set. (b) Model resolution
section for optimized data set generated by the ‘Compare R’ method. ............................ 84
Figure 62. Test inversion model for the different arrays. (a) Apparent resistivity
pseudosection for the simple dipole-dipole array with a dipole length ‘a’ of 1.0 meter
with the dipole separation factor ‘n’ of 1 to 6. (b) Apparent resistivity pseudosection for
the dipole-dipole array with overlapping data levels. (c) The synthetic model with 4
rectangular blocks of 100 .m embedded in a medium of 10 .m. as used by Wilkinson
et al. (2006b). .................................................................................................................... 85
Figure 63. Inversion results with the (a) Wenner array, (b) Wenner-Schlumberger array, (c)
simple dipole-dipole array and (d) dipole-dipole array with overlapping data levels. ..... 86
Figure 64. Inversions models with the (a) optimized data set with 4462 data points and (b) a
truncated optimized data set with 413 data points. ........................................................... 86
Figure 65. Comparison between the (a) model resolution, (b) model resolution index and (c)
DOI index sections for the Landfill data set. .................................................................... 87
Figure 66. The model resolution plots for the three streamer configurations. ......................... 88
Figure 67. An example of 3-D effects on a 2-D survey. (a) Apparent resistivity
pseudosections (Wenner array) along lines at different y-locations over (b) a 3-D
structure shown in the form of horizontal slices. .............................................................. 91
Figure 68. The 2-D sensitivity sections for the pole-dipole array with a dipole length of 1
meter and with (a) n=6, (b) n=12 and (c) n=18. Note that as the ‘n’ factor increases, the
zone of high positive sensitive values becomes increasingly concentrated in a shallower
zone below the P1-P2 dipole. ............................................................................................ 92
Figure 69. Example of apparent resistivity pseudosection with pole-dipole array with large
‘n’ values. Note that the anomaly due to a small near-surface high resistivity block
becomes greater as the ‘n’ factor increases. This means that the sensitivity of the array to
the near-surface region between the P1-P2 potential dipole becomes greater as the ‘n’
factor increases. ................................................................................................................. 93
Figure 70. Diagrammatic illustration of differences in objective function shapes for the pole-
pole array and dipole-dipole array data sets leading to different models obtained from
optimization routine. ......................................................................................................... 95
Figure 71. The I.P. values for some rocks and minerals. ......................................................... 97
Figure 72. The Cole-Cole model. (a) Simplified electrical analogue circuit model (after
Pelton et al. 1978). = resistivity, m = chargeability, = time constant, c = relaxation
constant. (b) Amplitude and phase response to sine wave excitation (frequency domain).
(c) Transient response to square wave current pulse (time domain). Most I.P. receivers
measure the integral of the decay voltage signal over a fixed interval, mt, as a measure of
the I.P. effect. .................................................................................................................... 97
Figure 73. Magusi River massive sulphide ore body inversion models. (a) Apparent resistivity
pseudosection. Resistivity inversion models obtained using the (b) perturbation and (c)
complex resistivity methods. (d) Apparent I.P. (metal factor) pseudosection. I.P. models
obtained using the (e) perturbation and (f) complex resistivity methods. ....................... 101
Figure 74. Sketch of separated cable spreads setup used (after Dahlin and Loke, 2015). ..... 102
Figure 75. Resistivity and chargeability pseudosection from field demo at 3rd IP workshop at
Ile d’Oleron (after Dahlin and Loke, 2015). .................................................................. 102
Figure 76. The possible arrangements of the electrodes for the pole-pole array in the cross-
borehole survey and the 2-D sensitivity sections. The locations of the two boreholes are
shown by the vertical black lines. ................................................................................... 104
Figure 77. A schematic diagram of two electrodes below the surface. The potential measured
at P can be considered as the sum of the contribution from the current source C and its
image C’ above the ground surface. ................................................................................ 105
Figure 78. The 2-D sensitivity patterns for various arrangements with the pole-bipole array.
The arrangement with (a) C1 and P1 in first borehole and P2 in second borehole, (b) C1
in the first borehole and both P1 and P2 in the second borehole, (c) all three electrodes in
the first borehole and (d) the current electrode on the ground surface. .......................... 106
Figure 79. The 2-D sensitivity patterns for various arrangements of the bipole-bipole array.
(a) C1 and P1 are in the first borehole, and C2 and P2 are in second borehole. (b) C1 and
C2 are in the first borehole, and P1 and P2 are in second borehole. In both cases, the
distance between the electrodes in the same borehole is equal to the separation between
the boreholes. The arrangements in (c) and (d) are similar to (a) and (b) except that the
distance between the electrodes in the same borehole is half the spacing between the
boreholes. ........................................................................................................................ 107
Figure 80. Possible measurement sequences using the bipole-bipole array. Other possible
measurements sequences are described in the paper by Zhou and Greenhalgh (1997). . 109
Figure 81. Several possible bipole-bipole configurations with a single borehole. (a) The C1
and C2 electrodes at depths of 3 and 2 meters respectively below the 0 meter mark. (b)
The C1 and P1 electrodes at depths of 3 and 2 meters respectively below the 0 meter
mark. (c). The C1 electrode is at a depth of 3 meters below the 0 meter mark while the
C2 electrode is on the surface. (d). The C1 electrode is at a depth of 3 meters below the 0
meter mark while the P1 electrode is on the surface. ...................................................... 111
Figure 82. A pole-bipole survey with a single borehole. The C1 electrode is at a depth of 3
meters below the 0 meter mark. ...................................................................................... 111
Figure 83. Test of optimized cross-borehole arrays with a synthetic model. (a) Two-layer test
model with conductive and resistive anomalies. Inversion models for (b) optimized data
set with all arrays, (c) 'standard' data set and (d) the reduced optimized data set that
excludes arrays with both current (or potential) electrodes in the same borehole. All the
data sets have 1875 data points. The outlines of the rectangular blocks showing their true
positions are also shown. ................................................................................................. 113
Figure 84. Schematic diagram of the MERIT method with the electrodes are planted along
the surface and directly below using the direct push technology (after Harro and Kruse,
2013)................................................................................................................................ 114
Figure 85. Inversion models for the different data sets for the data collected with electrodes at
surface and 7.62 m depth with 4 m horizontal spacing. Models for the (a) standard arrays
(405 data points), optimized arrays with (b) 403 and (c) 514 data points....................... 115
Figure 86. Landslide field example, Malaysia. (a) The apparent resistivity pseudosection for a
survey across a landslide in Cangkat Jering and (b) the interpretation model for the
subsurface. ....................................................................................................................... 116
Figure 87. Industrial pollution example, U.K. (a) The apparent resistivity pseudosection from
a survey over a derelict industrial site, and the (b) computer model for the subsurface. 117
Figure 88. Mapping of holes in a clay layer, U.S.A. (a) Apparent resistivity pseudosection for
the survey to map holes in the lower clay layer. (b) Inversion model and (c) sensitivity
values of the model cells used by the inversion program. .............................................. 118
Figure 89. Water infiltration mapping, U.K. (a) The apparent resistivity and (b) inversion
model sections from the survey conducted at the beginning of the Birmingham
infiltration study. This shows the results from the initial data set that forms the base
model in the joint inversion with the later time data sets. As a comparison, the model
obtained from the inversion of the data set collected after 10 hours of irrigation is shown
in (c). ............................................................................................................................... 119
Figure 90. Time-lapse sections from the infiltration study. The sections show the change in
the subsurface resistivity values with time obtained from the inversion of the data sets
collected during the irrigation and recovery phases of the study. ................................... 120
Figure 91. Hoveringham pumping test, U.K. (a) The apparent resistivity pseudosection at the
beginning of the test. The inversion model sections at the (b) beginning and (v) after 220
minutes of pumping. ........................................................................................................ 121
Figure 92. Percentage relative change in the subsurface resistivity values for the
Hoveringham pumping test. To highlight the changes in the subsurface resistivity, the
changes in the model resistivity are shown. Note the increase in the model resistivity
below the borehole with time. ......................................................................................... 121
Figure 93. Use of Archie’s Law for the Hoveringham pumping test. Sections showing the
relative desaturation values obtained from the inversion models of the data sets collected
during the different stages of the Hoveringham pumping test. Archie’s Law probably
gives a lower limit for the actual change in the aquifer saturation. ................................ 122
Figure 94. Groundwater survey, Nigeria. (a). Apparent resistivity pseudosection. (b) The
inversion model with topography. Note the location of the borehole at the 175 meters
mark................................................................................................................................. 123
Figure 95. The inversion model after 4 iterations from an underwater riverbed survey by
Sage Engineering, Belgium. ............................................................................................ 124
Figure 96. Thames River (CT, USA) survey with floating electrodes. (a) The measured
apparent resistivity pseudosection. Inversion models obtained (b) without constraints on
the water layer, and (c) with a fixed water layer. ............................................................ 125
Figure 97. Survey to map oil sands, Alberta, Canada. (a) Location of major tar sands deposits
in Alberta, Canada. (b) Example of resistivity log and geologic column of Athabasca oil
sands. (c) 2-D resistivity model from imaging survey (Kellett and Bauman, 1999). ..... 126
Figure 98. A simple arrangement of the electrodes for a 3-D survey. ................................... 128
Figure 99. Two possible measurement sequences for a 3-D survey. The location of potential
electrodes corresponding to a single current electrode in the arrangement used by (a) a
survey to measure the complete data set and (b) a cross-diagonal survey. ..................... 128
Figure 100. 3-D sensitivity plots for the pole-pole array. The plots are in the form of
horizontal slices through the earth at different depths. ................................................... 129
Figure 101. 3-D sensitivity plots for the pole-dipole array with n=1 in the form of horizontal
slices through the earth at different depths. The C1 electrode is the leftmost white cross.
......................................................................................................................................... 131
Figure 102. 3-D sensitivity plots for the pole-dipole array with n=4 in the form of horizontal
slices through the earth at different depths...................................................................... 131
Figure 103. 3-D sensitivity plots for the dipole-dipole array with n=1 in the form of
horizontal slices through the earth at different depths. The C2 electrode is the leftmost
white cross. ...................................................................................................................... 133
Figure 104. 3-D sensitivity plots for the dipole-dipole array with n=4 in the form of
horizontal slices through the earth at different depths. The C2 electrode is the leftmost
white cross. ...................................................................................................................... 133
Figure 105. The 3-D sensitivity plots for the Wenner alpha array at different depths. ......... 134
Figure 106. The 3-D sensitivity plots for the Wenner-Schlumberger array with the n=4 at
different depths................................................................................................................ 134
Figure 107. The 3-D sensitivity plots for the Wenner gamma array at different depths. ...... 135
Figure 108. Using the roll-along method to survey a 10 by 10 grid with a multi-electrode
system with 50 nodes. (a) Surveys using a 10 by 5 grid with the lines orientated in the x-
direction. (b) Surveys with the lines orientated in the y-direction. ................................. 136
Figure 109. A 3-D model with 4 rectangular prisms in a 15 by 15 survey grid. (a) The finite-
difference grid. (b) Horizontal apparent resistivity pseudosections for the pole-pole array
with the electrodes aligned in the x-direction.................................................................. 138
Figure 110. The model used in 3-D inversion. ...................................................................... 140
Figure 111. Inversion models for the Vetlanda landfill survey data set. (a) Using standard
inversion settings. (b) With a higher damping factor for the topmost layer. (c) Using
diagonal roughness filter in the horizontal (x-y) directions. (d) Using diagonal roughness
filters in the vertical (x-z and y-z) directions as well. (e) Using the roughness filter in all
directions. ........................................................................................................................ 142
Figure 112. Types of 3-D roughness filters. (a) With components in the x- and y- directions
only for the horizontal filter. (b) With components in the diagonal directions in the x-y
plane for the horizontal filter. (c) Applying the roughness filter with the corner model
cells as well. Only two (out of eight) corner cells are shown. ........................................ 143
Figure 113. Inversion models for the Vetlanda landfill survey data set using model cells of
equal lengths in the x- and y- directions. (a) Using standard inversion settings. (b) With a
higher damping factor for the topmost layer. (c) Using diagonal roughness filter in the
horizontal (x-y) directions. (d) Using diagonal roughness filters in the vertical (x-z and y-
z) directions as well. ........................................................................................................ 144
Figure 114. Synthetic model for long electrodes survey. 3-D model using cells of low
resistivity (0.01 .m) that are marked in red to simulate cased wells. ........................... 145
Figure 115. Comparison of inversion models using point electrodes with and without long
electrodes. (a) Inversion model for pole-pole data set using only surface point electrodes.
(b) Inversion model for pole-pole data set using 121 point and 3 long electrodes. ........ 146
Figure 116. 3-D data grid formats.......................................................................................... 147
Figure 117. Methods to model the effect of an electrode using the finite-difference and finite-
element methods.............................................................................................................. 148
Figure 118. The use of an appropriate mesh spacing to obtain sufficient accuracy for
electrodes in the same array that close together. ............................................................. 148
Figure 119. Map with survey lines and infrastructure at the Hanford site. ........................... 149
Figure 120. Types of model grids for the Hanford survey data set. (a) Using a 5 meters
spacing model for the entire area. (b) Using a 4 meters spacing model for the entire area.
(c) Using a mixture of 4 and 5 meters spacing model grid. ............................................ 150
Figure 121. Inversion model for Hanford survey site. ........................................................... 151
Figure 122. Arrangement of survey lines using a 3 cable system with the Abem SAS
instrument. ....................................................................................................................... 152
Figure 123. Inversions model for (a) combined Wenner-Schlumberger and dipole-dipole data
set, (b) optimized data set. The actual positions of the blocks are marked by black
rectangles......................................................................................................................... 152
Figure 124. Horizontal sections showing the model resolution for the (a) comprehensive data
set, (b) standard arrays and (c) optimized arrays. The electrode positions are marked by
small green crosses in the top layer in (a). ...................................................................... 153
Figure 125. Vertical cross-sections showing the model resolution for the comprehensive,
‘standard’, small and large optimized data sets............................................................... 154
Figure 126. The synthetic test model with two 1000 ohm.m rectangular blocks (marked by
black rectangles) embedded in a 100 ohm.m background medium. ............................... 154
Figure 127. The inversion models for the synthetic model with two 1000 ohm.m rectangular
blocks (marked by black rectangles) embedded in a 100 ohm.m background medium. 156
Figure 128. Comprehensive data set point-spread-function plots on a vertical x-z plane
located at y=4.5 m for a model cell at different depths with centre at x=4.5m and y=4.5m.
......................................................................................................................................... 156
Figure 129. Model resolution sections with circular perimeter for (a) comprehensive data set
with 180300 arrays, optimized data sets with (b) 946 and (c) 2000 arrays. .................... 157
Figure 130. Inversion results for survey with a circular perimeter using optimized arrays with
(a) 946 and (b) 2000 data points...................................................................................... 157
Figure 131. The (a) initial and (b) perturbed synthetic test models with apparent resistivity
pseudosections and inversion models assuming a constant electrode spacing and flat
surface. (c) The inversion model obtained with the algorithm that allows the electrodes
to shift.............................................................................................................................. 162
Figure 132. (a) 3-D synthetic test model with a rectangular survey grid. In the perturbed
model, the resistivities of the smaller blocks were changed from 400 and 20 ohm.m to
350 and 25 ohm.m. (b) The survey grid for the perturbed model. The four electrodes
shifted are marked by red circles. Electrode 1 was shifted 0.3 m in the x-direction,
electrode 2 moved 0.3 m in the y-direction, electrode 3 moved -0.2 m in both x and y-
directions while electrode 4 was shifted vertically upwards by 0.4 m. ........................... 163
Figure 133. Inversion models for the (a) initial and (b) perturbed data sets with fixed
electrodes in a rectangular grid. The true positions of the bands and prisms are marked by
black lines. The positions of the shifted electrodes are marked by small crosses in the top
layer in (b). ...................................................................................................................... 163
Figure 134. Inversion models for the perturbed data set using different relative damping
factors for the electrodes positions vector with a homogenous half-space starting model.
......................................................................................................................................... 164
Figure 135. Surface x-z profiles along the line y=10 m for inversions using a (a) homogenous
half-space and (b) initial data set starting models. .......................................................... 164
Figure 136. Inversion models for the perturbed data set using different relative damping
factors for the electrodes positions vector with the inversion model from the initial data
set as the starting model. ................................................................................................. 164
Figure 137. (a) Arrangement of electrodes in the Birmingham 3-D field survey. (b)
Horizontal and (c) vertical cross-sections of the model obtained from the inversion of the
Birmingham field survey data set. The locations of observed tree roots on the ground
surface are also shown..................................................................................................... 165
Figure 138. Example L-curve plots. (a) A plot of the model roughness versus the data misfit
for the Birmingham survey data set for a number of damping factor values. (b) A plot of
the curvature of the curve in (a) for the different damping factor values. The ‘optimum’
damping factor is at the maximum curvature point......................................................... 166
Figure 139. The 3-D model obtained from the inversion of the Lernacken Sludge deposit
survey data set displayed as horizontal slices through the earth. .................................... 167
Figure 140. A 3-D view of the model obtained from the inversion of the Lernacken Sludge
deposit survey data set displayed with the Slicer/Dicer program. A vertical exaggeration
factor of 2 is used in the display to highlight the sludge ponds. Note that the color contour
intervals are arranged in a logarithmic manner with respect to the resistivity. ............... 167
Figure 141. Map of the Panama Canal region with the survey area marked (Noonan and
Rucker, 2011). ................................................................................................................. 168
Figure 142. The model grid used for the Panama Canal floating electrodes survey with 20 by
20 meters cells. The location of the electrodes along the surveys lines are shown as
colored points. ................................................................................................................. 169
Figure 143. The inversion model for the Panama Canal floating electrodes survey. The first 4
layers correspond to the water column. ........................................................................... 169
Figure 144. Geological map of the Copper Hill area (White et al. 2001). ............................. 170
Figure 145. Electrodes layout used for the 3-D survey of the Copper Hill area.................... 170
Figure 146. The I.P. model obtained from the inversion of the Copper Hill survey data set.
Yellow areas have chargeability values of greater than 35 mV/V, while red areas have
chargeability values of greater than 45 mV/V (White et al., 2001). ............................... 171
Figure 147. Location of uranium mines in the Athabasca basin, Saskatchewan (Bingham et
al., 2006). ........................................................................................................................ 172
Figure 148. Geological model of uranium deposit (Bingham et al., 2006). .......................... 172
Figure 149. Example of inversion model from the Midwest deposit area showing the
resistivity at a depth of about 200 meters (Bingham et al., 2006). ................................. 173
Figure 150. A complex time-lapse field survey example. (a) Map of Cripple Creek survey
site. (b) Overhead view of the inversion model grid with electrodes layout. (c) Iso-surface
contours for the -4% resistivity change at different times after the injection of the sodium
cyanide solution (that started at 2.8 hours from the first data set in snapshots used). t1=
1.1 hours, t2= 2.4 hours, t3= 3.7 hours, t4= 4.9 hours. (d) Overhead view of iso-surfaces.
......................................................................................................................................... 174
Figure 151. Geological map of (a) south-west South Australia, (b) the Burra area, and (c) a
plot of survey electrodes and model cells layout. ........................................................... 175
Figure 152. Burra survey (a) resistivity and (b) I.P. inversion model layers. ........................ 176
Figure 153. The model (a) resistivity and (b) I.P. resolution index, and (c) VOI values. The
red arrows at the left side of (c) shows the position of the vertical slice shown in Figure
154. .................................................................................................................................. 177
Figure 154. Vertical cross-sections of the (a) resistivity model resolution index, (b) I.P
resolution index and (c) VOI in the X-Z plane 0.8 km north of the origin. .................... 178
Figure 155. 2-D test model with conductive overburden. ..................................................... 189
Figure 156. Models obtained with different relative damping weights αs. ............................ 190
Figure 157. Models for the Burra data set with a relative damping weight of 0.5 for αs. ..... 191
Figure 158. I.P. vertical sections along the y-direction (at x=1650 to 1700 m) Burra data set
(a) without (αs=0.0) and (b) with (αs=0.5) a reference model constraint. ....................... 191
Figure 159. A long electrode in a (a) homogeneous medium, (b) two-layer medium with a
low resistivity lower layer and (c) partly in air. .............................................................. 193
Figure 160. Example of non-conventional electrodes arrangements. (a) Offset pole-dipole
arrangement, (b) distributed pole-dipole arrangement. ................................................... 194
List of Tables
Table 1. 1-D inversion examples using the RES1D.EXE program. ........................................ 11
Table 2. The median depth of investigation (ze) for the different arrays (Edwards, 1977). L is
the total length of the array. Note identical values of ze/a for the Wenner-Schlumberger
and pole-dipole arrays. Please refer to Figure 4 for the arrangement of the electrodes for
the different arrays. The geometric factor is for an "a" value of 1.0 meter. For the pole-
dipole array, the array length ‘L’ only takes into account the active electrodes C1, P1 and
P2 (i.e. it does not take into account the remote C2 electrode). ........................................ 29
Table 3. Forward modeling examples. ..................................................................................... 50
Table 4. Methods to remove bad data points ........................................................................... 54
Table 5. Tests with different inversion options ....................................................................... 63
Table 6. Tests with different topographic modeling options ................................................... 68
Table 7. Tests with option to fix the model resistivity ............................................................ 70
Table 8. Inversion times for the long Redas survey data set using different settings. ............. 79
Table 9. Inversion times for different line lengths. .................................................................. 79
Table 10. Tests with 2-D I.P. inversion ................................................................................. 102
Table 11. A few borehole inversion tests............................................................................... 109
Table 12. 3-D forward modeling examples ........................................................................... 137
Table 13. 3-D inversion examples ......................................................................................... 141
1
where r is the distance of a point in the medium (including the ground surface) from the
electrode. In practice, all resistivity surveys use at least two current electrodes, a positive
current and a negative current source. Figure 2 show the potential distribution caused by a pair
of electrodes. The potential values have a symmetrical pattern about the vertical place at the
mid-point between the two electrodes. The potential value in the medium from such a pair is
given by
I 1 1
(1.7)
2 rC1 rC 2
where rC1 and rC2 are distances of the point from the first and second current electrodes.
In practically all surveys, the potential difference between two points (normally on the
ground surface) is measured. A typical arrangement with 4 electrodes is shown in Figure 3.
The potential difference is then given by
I 1 1 1 1 (1.8)
2 rC1P1 rC 2 P1 rC1P 2 rC 2 P 2
The above equation gives the potential that would be measured over a homogenous half space
with a 4 electrodes array.
Actual field surveys are conducted over an inhomogenous medium where the
subsurface resistivity has a 3-D distribution. The resistivity measurements are still made by
injecting current into the ground through the two current electrodes (C1 and C2 in Figure 3),
and measuring the resulting voltage difference at two potential electrodes (P1 and P2). From
the current (I) and potential ( ) values, an apparent resistivity ( a) value is calculated.
a k (1.9)
I
where 2
k
1 1 1 1
rC1P1 rC 2 P1 rC1P 2 rC 2 P 2
k is a geometric factor that depends on the arrangement of the four electrodes. Resistivity
measuring instruments normally give a resistance value, R = /I, so in practice the apparent
resistivity value is calculated by
a = k R (1.10)
The calculated resistivity value is not the true resistivity of the subsurface, but an
“apparent” value that is the resistivity of a homogeneous ground that will give the same
resistance value for the same electrode arrangement. The relationship between the “apparent”
resistivity and the “true” resistivity is a complex relationship. To determine the true subsurface
resistivity from the apparent resistivity values is the “inversion” problem. Methods to carry out
such an inversion will be discussed in more detail at the end of this chapter.
Figure 4 shows the common arrays used in resistivity surveys together with their
geometric factors. In a later section, we will examine the advantages and disadvantages of some
of these arrays.
There are two more electrical based methods that are closely related to the resistivity
method. They are the Induced Polarization (IP) method, and the Spectral Induced Polarization
(SIP) (also known as Complex Resistivity (CR)) method. Both methods require measuring
instruments that are more sensitive than the normal resistivity method, and with significantly
higher currents. IP surveys are comparatively more common, particularly in mineral
exploration surveys. It is able to detect conductive minerals of very low concentrations that
Figure 1. The flow of current from a point current source and the resulting potential
distribution.
Figure 2. The potential distribution caused by a pair of current electrodes. The electrodes are
1 m apart with a current of 1 ampere and a homogeneous half-space with resistivity of 1 m.
Figure 3. A conventional array with four electrodes to measure the subsurface resistivity.
Figure 4. Common arrays used in resistivity surveys and their geometric factors. Note that the
dipole-dipole, pole-dipole and Wenner-Schlumberger arrays have two parameters, the dipole
length “a” and the dipole separation factor “n”. While the “n” factor is commonly an integer
value, non-integer values can also be used.
Metals, such as iron, have extremely low resistivity values. Chemicals that are strong
electrolytes, such as potassium chloride and sodium chloride, can greatly reduce the resistivity
of ground water to less than 1 m even at fairly low concentrations. The effect of weak
electrolytes, such as acetic acid, is comparatively smaller. Hydrocarbons, such as xylene
(6.998x1016 m), typically have very high resistivity values. However, in practice the
percentage of hydrocarbons in a rock or soil is usually quite small, and might not have a
significant effect on the bulk resistivity. However when the concentration of the hydrocarbon
is high, such as the commercial oil sands deposits in Canada, the resistivity method has proved
to be a useful exploration method for such deposits (see section 7.10 for an example).
1.3 1-D resistivity surveys and inversions – applications, limitations and pitfalls
The resistivity method has its origin in the 1920’s due to the work of the Schlumberger
brothers. For approximately the next 60 years, for quantitative interpretation, conventional
sounding surveys (Koefoed, 1979) were normally used. In this method, the center point of the
electrode array remains fixed, but the spacing between the electrodes is increased to obtain
more information about the deeper sections of the subsurface.
The measured apparent resistivity values are normally plotted on a log-log graph paper.
To interpret the data from such a survey, it is normally assumed that the subsurface consists of
horizontal layers. In this case, the subsurface resistivity changes only with depth, but does not
change in the horizontal direction. A one-dimensional model of the subsurface is used to
interpret the measurements (Figure 6a). Figure 7 shows an example of the data from a sounding
survey and a possible interpretation model. This method has given useful results for geological
situations (such the water-table) where the one-dimensional model is approximately true.
The software provided, RES1D.EXE, is a simple inversion and forward modeling
program for 1-D models that consists of horizontal layers. In the software package, several files
with extensions of DAT are example data files with resistivity sounding data. Files with the
MOD extension are model files that can be used to generate synthetic data for the inversion
part of the program. As a first try, read in the file WENNER3.DAT that contains the Wenner
array sounding data for a simple 3-layer model.
Figure 6. The three different models used in the interpretation of resistivity measurements.
Figure 7. A typical 1-D model used in the interpretation of resistivity sounding data for the
Wenner array.
The greatest limitation of the resistivity sounding method is that it does not take into
account lateral changes in the layer resistivity. Such changes are probably the rule rather than
the exception. The failure to include the effect of such lateral changes can results in errors in
the interpreted layer resistivity and/or thickness. As an example, Figure 8 shows a 2-D model
where the main structure is a two-layer model with a resistivity of 10 m and a thickness of 5
meters for the upper layer, while the lower layer has a resistivity of 100 m. To the left of the
center point of the survey line, a low resistivity prism of 1 m is added in the upper layer to
simulate a lateral inhomogeneity. The 2-D model has 144 electrodes that are 1 meter apart. The
apparent resistivity pseudosections for the Wenner and Schlumberger array are also shown. For
the Schlumberger array, the spacing between the potential electrodes is fixed at 1.0 meter for
the apparent resistivity values shown in the pseudosection. The sounding curves that are
obtained with conventional Wenner and Schlumberger array sounding surveys with the mid-
point at the center of the line are shown in Figure 9. In the 2-D model, the low resistivity
rectangular prism extends from 5.5 to 18.5 meters to the left of the sounding mid-point. The
ideal sounding curves for both arrays for a two-layer model (i.e. without the low resistivity
prism) are also shown for comparison. For the Wenner array, the low resistivity prism causes
the apparent resistivity values in the sounding curve (Figure 9a) to be too low for spacing
values of 2 to 9 meters and for spacings larger than 15 meters. At spacings of between 9 to 15
meters, the second potential electrode P2 crosses over the low resistivity prism. This causes the
apparent resistivity values to approach the two-layer model sounding curve. If the apparent
resistivity values from this model are interpreted using a conventional 1-D model, the resulting
model could be misleading. In this case, the sounding data will most likely to be interpreted as
a three-layer model.
The effect of the low resistivity prism on the Schlumberger array sounding curve is
slightly different. The apparent resistivity values measured with a spacing of 1 meter between
the central potential electrodes are shown by black crosses in Figure 9b. For electrode spacings
(which is defined as half the total length of the array for the Schlumberger array) of less than
15 meters, the apparent resistivity values are less than that of the two-layer sounding curve. For
spacings greater than 17 meters, the apparent resistivity values tend to be too high. This is
probably because the low resistivity prism lies to the right of the C2 electrode (i.e. outside the
array) for spacings of less than 15 meters. For spacings of greater than 17 meters, it lies between
the P2 and C2 electrodes. Again, if the data is interpreted using a 1-D model, the results could
be misleading. One method that has been frequently recommended to “remove” the effect of
lateral variations with the Schlumberger array is by shifting curve segments measured with
different spacings between the central potential electrodes. The apparent resistivity values
measured with a spacing of 3 meters between the potential electrodes are also shown in Figure
9b. The difference in the sounding curves with the spacings of 1 meter and 3 meters between
the potential electrodes is small, particularly for large electrode spacings. Thus any shifting in
the curve segments would not remove the distortion in the sounding curve due to the low
resistivity prism. The method of shifting the curve segments is probably more applicable if the
inhomogeneity lies between the central potential electrodes, and probably ineffective if the
inhomogeneity is beyond the largest potential electrodes spacing used (which is the case in
Figure 8). However, note that the effect of the prism on the Schlumberger array sounding curve
is smaller at the larger electrode spacings compared with the Wenner array (Figure 9). The
main reason is probably the larger distance between the P2 and C2 electrodes in the
Schlumberger array.
A more reliable method to reduce the effect of lateral variations on the sounding data
is the offset Wenner method (Barker, 1978). It makes use of the property that the effect of an
inhomogeneity on the apparent resistivity value is of opposite sign if it lies between the two
potential electrodes or if it is between a potential and a current electrode. For the example
shown in Figure 8, if the low resistivity body lies in between a current and potential electrode
(the P2 and C2 electrodes in this case), the measured apparent resistivity value would be lower.
If the low resistivity body lies in between the P1 and P2 electrodes, it will cause the apparent
resistivity value to be higher. The reason for this phenomenon can be found in the sensitivity
pattern for the Wenner array (see Figure 23a). By taking measurements with different positions
for the mid-point of the array, the effect of the low resistivity body can be reduced.
Figure 8. A 2-D two-layer model with a low resistivity prism in the upper layer. The calculated
apparent resistivity pseudosections for the (a) Wenner and (b) Schlumberger arrays. (c) The
2D model. The mid-point for a conventional sounding survey is also shown.
Another classical survey technique is the profiling method. In this case, the spacing
between the electrodes remains fixed, but the entire array is moved along a straight line. This
gives some information about lateral changes in the subsurface resistivity, but it cannot detect
vertical changes in the resistivity. Interpretation of data from profiling surveys is mainly
qualitative.
The most severe limitation of the resistivity sounding method is that horizontal (or
lateral) changes in the subsurface resistivity are commonly found. The ideal situation shown in
Figure 6a is rarely found in practice. As shown by the examples in Figure 8 and Figure 9,
lateral changes in the subsurface resistivity will cause changes in the apparent resistivity values
that might be, and frequently are, misinterpreted as changes with depth in the subsurface
resistivity. In many engineering and environmental studies, the subsurface geology is very
complex where the resistivity can change rapidly over short distances. The 1-D resistivity
Figure 9. Apparent resistivity sounding curves for a 2-D model with a lateral inhomogeneity.
(a) The apparent resistivity curve extracted from the 2D pseudosection for the Wenner array.
The sounding curve for a two-layer model without the low resistivity prism is also shown by
the black line curve. (b) The apparent resistivity curves extracted from the 2-D pseudosection
for the Schlumberger array with a spacing of 1.0 meter (black crosses) and 3.0 meters (red
crosses) between the potential electrodes. The sounding curve for a two-layer model without
the low resistivity prism is also shown.
To use the RES1D.EXE program for the exercises in the table below, as well as the
other programs that we shall use in the later sections, follows the usual sequence used by
Windows XP/Vista/7/8/10. Click the ‘Start’ button, followed by ‘Programs’ and the look for
the RES1D folder in the list of installed programs. Alternatively, you can create a shortcut icon
on the Windows Desktop.
To obtain a more accurate subsurface model than is possible with a simple 1-D model,
a more complex model must be used. In a 2-D model (Figure 6b), the resistivity values are
allowed to vary in one horizontal direction (usually referred to as the x direction) but assumed
to be constant in the other horizontal (the y) direction. This approximation is reasonable for
survey lines that are perpendicular to the strike of an elongated structure. The most realistic
model would be a fully 3-D model (Figure 6c) where the resistivity values are allowed to
change in all 3 directions. The use of 2-D and 3-D surveys and interpretation techniques will
be examined in detail in the following chapters.
where m is the number of measurements. The model response f can be written in a similar form.
f col( f1 , f 2 ,....., f m ) (1.13)
For resistivity problems, it is a common practice to use the logarithm of the apparent resistivity
values for the observed data and model response, and the logarithm of the model values as the
model parameters. The model parameters can be represented by the following vector
q col(q1 , q2 ,....., qn ) (1.14)
where n is the number of model parameters. The difference between the observed data and the
model response is given by the discrepancy vector g that is defined by
g=y-f (1.15)
In the least-squares optimization method, the initial model is modified such that the
sum of squares error E of the difference between the model response and the observed data
values is minimized.
n
E g Tg g i2 (1.16)
i 1
To reduce the above error value, the following Gauss-Newton equation is used to
determine the change in the model parameters that should reduce the sum of squares error
(Lines and Treitel 1984).
J T J Δqi J T g (1.17)
where q is the model parameter change vector, and J is the Jacobian matrix (of size m by n)
of partial derivatives. The elements of the Jacobian matrix are given by
f
J ij i (1.18)
q j
that is the change in the ith model response due to a change in the jth model parameter. After
calculating the parameter change vector, a new model is obtained by
q k 1 q k Δqk (1.19)
In practice, the simple least-squares equation (1.17) is rarely used by itself in
geophysical inversion. In some situations the matrix product J T J might be singular, and thus
the least-squares equation does not have a solution for q. Another common problem is that
the matrix product J T J is nearly singular. This can occur if a poor initial model that is very
different from the optimum model is used. The parameter change vector calculated using
equation (1.17) can have components that are too large such that the new model calculated with
(1.19) might have values that are not realistic. One common method to avoid this problem is
the Marquardt-Levenberg modification (Lines and Treitel, 1984) to the Gauss-Newton
equation that is given by
J T
J I Δqk J T g (1.20)
where I is the identity matrix. The factor is known as the Marquardt or damping factor, and
this method is also known as the ridge regression method (Inman 1975) or damped least-
squares method. The damping factor effectively constrains the range of values that the
components of parameter change vector can q take. While the Gauss-Newton method in
equation (1.17) attempts to minimize the sum of squares of the discrepancy vector only, the
Marquardt-Levenberg method also minimizes a combination of the magnitude of the
discrepancy vector and the parameter change vector. This method has been successfully used
in the inversion of resistivity sounding data where the model consists of a small number of
layers. For example, it was used in the inversion of the resistivity sounding example in Figure
7 with three layers (i.e. five model parameters). However when the number of model
parameters is large, such as in 2-D and 3-D inversion models that consist of a large number of
small cells, the model produced by this method can have an erratic resistivity distribution with
spurious high or low resistivity zones (Constable et al., 1987). To overcome this problem, the
Gauss-Newton least-squares equation is further modified so as to minimize the spatial
variations in the model parameters (i.e. the model resistivity values change in a smooth or
gradual manner). This smoothness-constrained least-squares method (Ellis and Oldenburg
1994a, Loke 2011) has the following mathematical form.
J T J F Δqk J T g Fq k , (1.21)
where F x C x C x y C y C y z C z C z
T T T
and Cx, Cy and Cz are the roughness filter matrices in the x-, y- and z-directions that couples
the model blocks in those directions (Figure 10, Figure 112). x, y and z are the relative
weights given to the roughness filters in the x-, y- and z-directions. One common form of the
roughness filter matrix is the first-order difference matrix (deGroot-Hedlin and Constable
1990) that is given by
1 1 0 0 .. .. .. 0
0 1 1 0 .. .. .. 0
0 0 1 1 0 .. .. 0
..
C (1.22)
..
..
..
0
Figure 10. Coupling between neighboring model cells through the roughness filter in a 2-D
model. (a) In the horizontal and vertical diretcions only, and (b) in diagonal diretcions as well.
Equation 1.21 also tries to minimize the square of the spatial changes, or roughness, of
the model resistivity values. It is in fact an l2 norm smoothness-constrained optimization
method. This tends to produce a model with a smooth variation of resistivity values. This
approach is acceptable if the actual subsurface resistivity varies in a smooth or gradational
manner. In some cases, the subsurface geology consists of a number of regions that are
internally almost homogeneous but with sharp boundaries between different regions. For such
cases, the inversion formulation in (1.21) can be modified so that it minimizes the absolute
changes in the model resistivity values (Claerbout and Muir, 1973). This can sometimes give
significantly better results. Technically this is referred to as an l1 norm smoothness-constrained
optimization method, or more commonly known as the blocky inversion method. A number of
techniques can be used for such a modification. One simple method to implement an l1 norm
based optimization method using the standard least-squares formulation is the iteratively
reweighted least-squares method (Wolke and Schwetlick, 1988; Farquharson and Oldenburg,
1998). The optimization equation in (1.21) is modified to
J T
R d J FR Δqk J T R d g FR q k , (1.23)
with FR x C Tx R m C x y C Ty R m C y z C Tz R m C z
where Rd and Rm are weighting matrices introduced so that different elements of the data misfit
and model roughness vectors are given equal weights in the inversion process.
Equation (1.23) provides a general method that can be further modified if necessary to
include known information about the subsurface geology. As an example, if it is known that
the variations in the subsurface resistivity are likely to be confined to a limited zone, the
damping factor values can modified (Ellis and Oldenburg 1994a) such that greater changes are
allowed in that zone. If the errors in the data points are known, a diagonal weighting matrix
can be used to give greater weights to data points with smaller errors.
One important modification to the least-squares optimization method is in time-lapse
inversion using the following equation (Kim et al. 2009, Loke et al. 2014a).
J T
i
R d J i λ i FR α M T R t M Δri J Ti R d g i λ i FR α i M T R t M r i 1 (1.24)
M is the difference matrix applied across the time models with only the diagonal and one sub
diagonal elements having values of 1 and -1, respectively, while r is the combined resistivity
model for all the time series. is the temporal damping factor that gives the weight for
minimizing the temporal changes in the resistivity compared to the model roughness and data
misfit. Note the smoothness-constraint is not only applied in space through the F matrix, but
also across the different time models through the M matrix. This is illustrated in Figure 11.
One common variation is to apply a constraint such that the model is ‘close’ to a
reference model, qR. Equation (1.24) then then modified to the following form.
J T J FR Δqk J T R dg (FR I) q k q R (1.25)
qR is usually a homogeneous background model. This equation imposes and additional
constraint the new model is close to the reference model with the damping factor weight α. Its
effect is similar to the Marquardt constraint in equation 1.20, and it prevents very large
deviations from the reference model.
The smoothness-constrained least-squares optimization method involves the damping
factor term . This term balances the need to reduce the data misfit (so that the calculated
apparent resistivity values are as close as possible to the measured values) while producing a
model that is ‘reasonably’ smooth and perhaps more geologically realistic. A smaller value of
term will generally produce a model with a lower data misfit but usually at the expense of
larger variations in the model resistivity values. If the error of the measured apparent resistivity
values is known, a prudent approach might be to select a model where the data misfit is similar
to the known measurement errors. However, for most field data sets, the measurement error is
not known. There are generally two methods to automatically select the ‘optimum’ damping
factor for such cases, the GCV and L-curve methods (Farquharson and Oldenburg, 2004). For
the inversion of a single data set, either method can be used. However, for time-lapse data sets
(Loke et al. 2014a), the L-curve method has a clear advantage due to the sparse block structure
of the Jacobian matrix (Figure 11b). The GCV method requires a matrix inversion. However,
the inverse of a sparse matrix is usually a full matrix. This makes it impractical to use for time-
lapse models with many time-series measurements such that the combined model might have
hundreds of thousands of model parameters (i.e. the number of model parameters for a single
data set multiplied by the number of time series data sets).
Figure 11. (a) Coupling between corresponding model blocks in two time-lapse models using
a cross-model time-lapse smoothness constraint. (b) Example Jacobian matrix structure for five
time series data sets and models. Each grey rectangle represents the Jacobian matrix associated
with a single set of measurements.
lateral changes, they are often mistakenly modeled as changes in the depths of the boundaries.
More recent efforts have been in combining the cell based and boundary based
inversion methods (Smith et al., 1999). One such method is the laterally constrained inversion
method (Auken and Christiansen, 2004). In this method, lateral changes (but not vertical
changes) are allowed in each region (Figure 12c), and abrupt transitions across the boundaries
are also allowed. The “model parameters” for the example in Figure 12 are then the twenty-
four resistivity values (1 and 24) and the depths at thirteen points (z1 to z13) along the boundary
giving a total of thirty-seven parameters. Information from other sources, such as borehole or
seismic data, can be used to provide an initial estimate of the depth to the boundary. A common
situation is when the depth information is available at only one borehole. In this case, the initial
boundary is usually set a constant depth. The inversion method then adjusts the depths at a
number of points along the boundary during the inversion process. A smoothness-constraint is
applied to minimize changes in the depths between adjacent points on the same boundary
(Smith et al., 1999). This method works particularly well where the subsurface consists of
several sedimentary zones.
A further generalization of this concept is to allow both vertical and lateral changes
within each region (as in a pure cell based model) while also allowing sharp changes across the
boundaries (Figure 12d). The model shown in Figure 12d has seventy-two resistivity values
and five depth values, giving a total of seventy-seven model parameters. This type of
discretization is particularly useful where near-surface inhomogeneities that occur at different
depths within the top layer have a large effect on the measured apparent resistivity values
(Figure 13). The model cells used so far has rectangular shapes (Figure 12). This is partly due
to the use of the finite-difference method in calculating the model apparent resistivity values.
A slight disadvantage is that the boundary is approximated by a series of rectangular steps.
Figure 12e shows a possible variation using the finite-element method with trapezoidal cells
where the edges of the cells adjacent to the boundary are adjusted so as to conform to the true
shape of the boundary.
Figure 12. The different models for the subsurface used in the interpretation of data from 2-D
electrical imaging surveys. (a) A purely cell based model. (b) A purely boundary based model.
(c) The laterally constrained model. (d) A combined cell based and boundary based model with
rectangular cells, and (e) with boundary conforming trapezoidal cells.
Figure 13. Example of a cell based model with a variable boundary. (a) The test model. (b)
The apparent resistivity pseudosection. (c) The inversion model with a variable sharp boundary
that is marked by a black line.
2.1 Introduction
We have seen that the greatest limitation of the resistivity sounding method is that it
does not take into account horizontal changes in the subsurface resistivity. A more accurate
model of the subsurface is a two-dimensional (2-D) model where the resistivity changes in the
vertical direction, as well as in the horizontal direction along the survey line. It is assumed that
resistivity does not change in the direction that is perpendicular to the survey line. In many
situations, particularly for surveys over elongated geological bodies, this is a reasonable
assumption. In theory, a 3-D resistivity survey and interpretation model should be even more
accurate. However, at the present time, 2-D surveys are the most practical economic
compromise between obtaining very accurate results and keeping the survey costs down
(Dahlin, 1996). Typical 1-D resistivity sounding surveys usually involve about 10 to 20
readings, while 2-D imaging surveys involve about 100 to 1000 measurements. In comparison,
a 3-D survey might involve several thousand measurements.
The cost of a typical 2-D survey could be several times the cost of a 1-D sounding
survey, and is probably comparable with a seismic refraction survey. In many geological
situations, 2-D electrical imaging surveys can give useful results that are complementary to the
information obtained by other geophysical methods. For example, seismic methods can map
undulating interfaces well, but will have difficulty (without using advanced data processing
techniques) in mapping discrete bodies such as boulders, cavities and pollution plumes. Ground
radar surveys can provide more detailed pictures but have very limited depth penetration in
areas with conductive unconsolidated sediments, such as clayey soils. Two-dimensional
electrical surveys can be used in conjunction with seismic or GPR surveys as they provide
complementary information about the subsurface.
Wenner electrode array for a system with 20 electrodes. In this example, the spacing between
adjacent electrodes is “a”. The first step is to make all the possible measurements with the
Wenner array with an electrode spacing of “1a”. For the first measurement, electrodes number
1, 2, 3 and 4 are used. Notice that electrode 1 is used as the first current electrode C1, electrode
2 as the first potential electrode P1, electrode 3 as the second potential electrode P2 and
electrode 4 as the second current electrode C2. For the second measurement, electrodes number
2, 3, 4 and 5 are used for C1, P1, P2 and C2 respectively. This is repeated down the line of
electrodes until electrodes 17, 18, 19 and 20 are used for the last measurement with “1a”
spacing. For a system with 20 electrodes, note that there are 17 (20 - 3) possible measurements
with “1a” spacing for the Wenner array.
After completing the sequence of measurements with “1a” spacing, the next sequence
of measurements with “2a” electrode spacing is made. First electrodes 1, 3, 5 and 7 are used
for the first measurement. The electrodes are chosen so that the spacing between adjacent
electrodes is “2a”. For the second measurement, electrodes 2, 4, 6 and 8 are used. This process
is repeated down the line until electrodes 14, 16, 18 and 20 are used for the last measurement
with spacing “2a”. For a system with 20 electrodes, note that there are 14 (20 - 2x3) possible
measurements with “2a” spacing.
Figure 14. The arrangement of electrodes for a 2-D electrical survey and the sequence of
measurements used to build up a pseudosection.
The same process is repeated for measurements with “3a”, “4a”, “5a” and “6a”
spacings. To get the best results, the measurements in a field survey should be carried out in a
systematic manner so that, as far as possible, all the possible measurements are made. This will
affect the quality of the interpretation model obtained from the inversion of the apparent
resistivity measurements (Dahlin and Loke, 1998).
Note that as the electrode spacing increases, the number of measurements decreases.
The number of measurements that can be obtained for each electrode spacing, for a given
number of electrodes along the survey line, depends on the type of array used. The Wenner
array gives the smallest number of possible measurements compared to the other common
arrays that are used in 2-D surveys.
The survey procedure with the pole-pole array is similar to that used for the Wenner
array. For a system with 20 electrodes, firstly 19 of measurements with a spacing of “1a” are
made, followed by 18 measurements with “2a” spacing, followed by 17 measurements with
“3a” spacing, and so on.
For the dipole-dipole, Wenner-Schlumberger and pole-dipole arrays (Figure 1.4), the
survey procedure is slightly different. As an example, for the dipole-dipole array, the
measurement usually starts with a spacing of “1a” between the C1-C2 (and also the P1-P2)
electrodes. The first sequence of measurements is made with a value of 1 for the “n” factor
(which is the ratio of the distance between the C1-P1 electrodes to the C1-C2 dipole length),
followed by “n” equals to 2 while keeping the C1-C2 dipole pair spacing fixed at “1a”. When
“n” is equals to 2, the distance of the C1 electrode from the P1 electrode is twice the C1-C2
dipole length. For subsequent measurements, the “n” spacing factor is usually increased to a
maximum value of about 6, after which accurate measurements of the potential are difficult
due to very low potential values. To increase the depth of investigation, the spacing between
the C1-C2 dipole pair is increased to “2a”, and another series of measurements with different
values of “n” is made. If necessary, this can be repeated with larger values of the spacing of
the C1-C2 (and P1-P2) dipole pairs. A similar survey technique can be used for the Wenner-
Schlumberger and pole-dipole arrays where different combinations of the “a” spacing and “n”
factor can be used.
It should be noted that in practice the sequence of measurements to take should be
arranged so that electrode polarization effects do not occur (Dahlin, 2000). This happens when
an electrode that was used as a current electrode is used as potential electrode in a following
measurement within a short time.
One technique used to extend horizontally the area covered by the survey, particularly
for a system with a limited number of electrodes, is the roll-along method. After completing
the sequence of measurements, the cable is moved past one end of the line by several unit
electrode spacings. All the measurements that involve the electrodes on part of the cable that
do not overlap the original end of the survey line are repeated (Figure 15).
Figure 15. The use of the roll-along method to extend the area covered by a 2-D survey.
Most of the above manufacturers have sub-agents in different countries, and there are
probably a few others that are not on the list. An interesting recent development is the
availability of lower cost Russian and Chinese systems. A few academic and research
institutions have designed their own systems, for example the British Geological Survey has
built an electrostatic based mobile system as well as an automatic system for time-lapse
surveys.
The instrument type can be divided into two broad categories, static and dynamic
systems. Most instruments are of the static type where many electrodes are connected to a
multi-electrode cable and planted into the ground during the survey. A typical static system is
the Abem Lund system shown in Figure 16. One common configuration is a split spread type
of cable connection to the switching unit at the center to reduce the individual cable length and
weight. The weight of a cable roll is directly proportional to the number of nodes and the
spacing between the nodes! A common spacing used for most engineering and environmental
surveys is 5 meters. Most systems come with a minimum of 28 nodes, with some system having
up to 128 nodes or more! The Lund system is a little unusual in that there are 4 individual
cables. Most systems use a 2 cables arrangement. The static systems can be further divided into
two sub-categories depending on the arrangement for the switching of the electrodes. Most of
the systems house the switches in a single unit and uses a cable with many individual wires
connected to each node. Typical examples are the Abem Lund and Campus systems. Another
arrangement is to have a small switching unit at each electrode and a cable with the minimum
number of wires. One early example is the Campus MRT system (Griffiths and Turnbull,
1985). More recent examples are the PASI and Iris Syscal systems.
There have been two new and interesting developments in the resistivity meter systems.
One is the addition of I.P. capability. However, the usefulness of the I.P. data from most multi-
electrode systems is of limited use. Most systems use a battery power source that cannot deliver
enough current (usually less than 1 ampere) for reliable I.P. signals, so the I.P. data is often
extremely noisy. However, there have been recent developments to reduce the noise in the I.P.
measurements (Dahlin and Leroux, 2012).
The second new development is multi-channel measuring systems. In such a system, a
number of potential measurements can be simultaneously made for a single pair of current
electrodes. This could significantly reduce the survey time. With certain array configurations,
a single 2-D survey line could involve thousands of measurements. The major part of the survey
time is waiting for a single channel instrument to complete the measurements that could take
more than several hours! The IP and multi-channel capability are relatively new developments,
so you will need to check the manufacturer's web site to get the latest information.
Figure 16. Sketch outline of the ABEM Lund Imaging System. Each mark on the cables
indicates an electrode position (Dahlin, 1996). The cables are placed along a single line (the
sideways shift in the figure is only for clarity). This figure also shows the principle of moving
cables when using the roll-along technique.
Static systems use a large number of nodes to get a wide data coverage. In contrast,
dynamic systems use a small number of nodes but move the entire system to obtain a wide
coverage. An example of such a system designed by Aarhus University in Denmark (Sorenson,
1996) is shown in Figure 17. A 100 meters cable with nine heavy cylindrical electrodes is
pulled by a small vehicle. Two of the electrodes are used as current electrodes, while six of
them are used for the potential measurements and one is used as a ground electrode. This
system relies on the current being injected into the ground by direct contact, so it can only be
used in open ground, such as farmlands in Northern Europe. A Wenner-Schlumberger type of
arrangement (Figure 4) is used but with non-integer “n” values for some of the measurements.
Another mobile system that does not require direct contact with the ground but uses capacitive
coupling (Gerard and Tabbagh, 1991; Shima et al., 1996, Panissod et al., 1998; Loke et al.,
2011a) to induce the flow of current in the ground. This system can be used in areas that are
paved, such as roads and city areas. One such system shown in Figure 18 is the Geometrics
OhmMapper system where a cable with 4 to 6 electrodes is attached to a measuring unit is
pulled by a single operator. The dipole-dipole type of arrangement is used but with non-integer
“n” values for some measurements.
One of main problems faced by mobile systems on land is to get sufficient current to
flow into the ground. Direct contact systems such as the Aarhus Pulled Array System can only
be used in areas with open ground. The capacitive coupling type does not require direct ground
contact and thus can be used in many areas where normal resistivity surveying systems cannot
be used (for example in built-up areas) but has the problem of a more limited depth of
penetration due to the limited amount of current that can be induced into the ground compared
to direct contact systems. An underwater environment provides an almost ideal situation for a
direct contact type of mobile system since there is no problem in obtaining good electrode
contact! Figure 19 shows a possible arrangement for an underwater mobile surveying system
where a cable with a number of nodes is pulled along the river/lake/sea bottom by a boat. Two
of the nodes are used as current electrodes, while the rest are used as potential electrodes. An
example of such an underwater survey is described in section 7.9. If this system is coupled with
a multi-channel resistivity meter, the survey can be carried out very rapidly. Shallow seismic
reflection surveys are frequently used in rivers/lakes/marine environments for engineering site
surveys. A mobile resistivity survey might be a useful addition in some situations, such as in
seismically opaque areas. In theory, both surveys can be carried out simultaneously to reduce
costs. The arrangement in Figure 19 has the cable dragged along the sea or river bed. This
places the electrodes close to the targets of interest but has the disadvantage that the cable can
get snagged be obstructions on the bottom. Another system uses a cable floating on the water
surface (section 7.9) but this has the disadvantage that most of the current flows within the
water layer, and very little goes into the sub-bottom materials which are being mapped. To
avoid both problems, some recent surveys use a submerged cable that is suspended about a
meter above the water bottom.
pulling vehicle
c c
c
current elec.
potential elec.
Figure 17. The Aarhus Pulled Array System. The system shown has two current (C) electrodes
and six potential electrodes (Christensen and Sørensen 1998, Bernstone and Dahlin 1999).
Figure 18. The Geometrics OhmMapper system using capacitive coupled electrodes.
(Courtesy of Geometrics Inc.).
Figure 19. Schematic diagram of a possible mobile underwater survey system. The cable has
two fixed current electrodes and a number of potential electrodes so that measurements can be
made at different spacings. The above arrangement uses the Wenner-Schlumberger type of
configuration. Other configurations, such as the gradient array, can also be used.
Figure 20. The apparent resistivity pseudosections from 2-D imaging surveys with different
arrays over a rectangular prism.
Figure 21. The parameters for the sensitivity function calculation at a point (x,y,z) within a half-
space. A pole-pole array with the current electrode at the origin and the potential electrode “a”
meters away is shown.
For the special case of a homogeneous half-space, the potential at a point in the half-
space due to a unit current source on the surface has a relatively simple form, which is
, and similarly
2 x 2 y 2 z 2
0.5
' (2.2)
2 x a y 2 z 2
2
0.5
After differentiating the above equations to obtain the divergence, and substituting into (2.1)
we get
1 x x a y 2 z 2
dxdydz (2.3)
V 4 2 x 2 y 2 z 2 1.5 x a 2 y 2 z 2
1.5
The 3-D Frechet derivative is then given by the term within the integral, i.e.
1 x x a y 2 z 2
F3D x, y, z (2.4)
4 2 x 2 y 2 z 2 1.5 x a 2 y 2 z 2
1.5
This gives the Frechet derivative or sensitivity function for the pole-pole array consisting of
just one current and one potential electrode. To obtain the Frechet derivative for a general four
electrodes array, we need to just add up the contributions from the four current-potential pairs,
just as we have done earlier for the potential in equation (1.8).
section of the earth above the "median depth of investigation" has the same influence on the
measured potential as the lower section. This tells us roughly how deep we can see with an
array. This depth does not depend on the measured apparent resistivity or the resistivity of the
homogeneous earth model. It should be noted that the depths are strictly only valid for a
homogeneous earth model, but they are probably good enough for planning field surveys. If
there are large resistivity contrasts near the surface, the actual depth of investigation could be
somewhat different.
The sensitivity function for other arrays can be determined by adding up the
contributions from the appropriate four pairs of current-potential electrodes. Figure 22b shows
the sensitivity function plot for the Wenner (alpha) array. Note that the curve around the
maximum is narrower for the Wenner array compared with the pole-pole array. This implies
that the Wenner array has a better vertical resolution than the pole-pole array.
Table 2 gives the median depth of investigation for the different arrays. To determine
the maximum depth mapped by a particular survey, multiply the maximum “a” electrode
spacing, or maximum array length “L“, by the appropriate depth factor given in Table 2. For
example, if the maximum electrode “a” spacing used by the Wenner array is 100 meters (or
maximum L 300 meters), then the maximum depth mapped is about 51 meters. For the dipole-
dipole, pole-dipole and Wenner-Schlumberger arrays, the “n” factor (Figure 4) must also be
taken into consideration. For the arrays with four active electrodes (such as the dipole-dipole,
Wenner and Wenner-Schlumberger arrays), it is probably easier to use the total array length
“L”. As an example, if a dipole-dipole survey uses a maximum value of 10 meters for “a” and
a corresponding maximum value of 6 for n, then the maximum “L” value is 80 meters. This
gives a maximum depth of investigation of 80x0.216 or about 17 meters.
Table 2 also includes the geometric factor for the various arrays for an "a" spacing of
1.0 meter. The inverse of the geometric factor gives an indication of the voltage that would be
measured between the P1 and P2 potential electrodes. The ratio of this potential compared to
the Wenner alpha array is also given, for example a value of 0.01 means that the potential is
1% of the potential measured by the Wenner alpha array with the same "a" spacing.
Figure 22. A plot of the 1-D sensitivity function. (a) The sensitivity function for the pole-
pole array. Note that the median depth of investigation (red arrow) is more than twice the
depth of maximum sensitivity (blue arrow). (b) The sensitivity function and median depth of
investigation for the Wenner array.
Table 2. The median depth of investigation (ze) for the different arrays (Edwards, 1977). L is
the total length of the array. Note identical values of ze/a for the Wenner-Schlumberger and
pole-dipole arrays. Please refer to Figure 4 for the arrangement of the electrodes for the
different arrays. The geometric factor is for an "a" value of 1.0 meter. For the pole-dipole array,
the array length ‘L’ only takes into account the active electrodes C1, P1 and P2 (i.e. it does not
take into account the remote C2 electrode).
z e/a z e/L Geometric Inverse Geometric
Factor Factor (Ratio)
Wenner Alpha 0.519 0.173 6.2832 0.15915 (1.0000)
Wenner Beta 0.416 0.139 18.850 0.05305 (0.3333)
Wenner Gamma 0.594 0.198 9.4248 0.10610 (0.6667)
Equatorial dipole-dipole
n=1 0.451 0.319 21.452 0.04662 (0.2929)
n=2 0.809 0.362 119.03 0.00840 (0.0528)
n=3 1.180 0.373 367.31 0.00272 (0.0171)
n=4 1.556 0.377 841.75 0.00119 (0.0075)
Wenner - Schlumberger
n=1 0.519 0.173 6.2832 0.15915 (1.0000)
n=2 0.925 0.186 18.850 0.05305 (0.3333)
n=3 1.318 0.189 37.699 0.02653 (0.1667)
n=4 1.706 0.190 62.832 0.01592 (0.1000)
n=5 2.093 0.190 94.248 0.01061 (0.0667)
n=6 2.478 0.191 131.95 0.00758 (0.0476)
n=7 2.863 0.191 175.93 0.00568 (0.0357)
n=8 3.247 0.191 226.19 0.00442 (0.0278)
n=9 3.632 0.191 282.74 0.00354 (0.0222)
n = 10 4.015 0.191 345.58 0.00289 (0.0182)
2.5.2b The median depth of investigation and the Gamma array problem
The use of the median depth of investigation runs into an interesting problem with the
Gamma array type (Carpenter and Habberjam, 1956) for certain configurations. It gives
ridiculously large and even infinite values for certain electrode arrangements (Barker, 1989).
This is because the sensitivity function changes from positive to negative with the depth.
k
2
2
0.5
for x>0.5a,
2 x 2 z 2 , 2 x a z 2 , xa
2
sensitivity values are shown from a depth of 0.025 meter downwards to 1.0 meter. In all the
sensitivity section diagrams, the location of the plotting point used in the pseudosection is
marked by a small black cross.
Figure 23. 2-D sensitivity sections for the Wenner array. The sensitivity sections for the (a)
alpha, (b) beta and (c) gamma configurations.
Figure 24 shows the sensitivity sections for this array for "n" values ranging from 1 to
6. The largest sensitivity values are generally located between the C2-C1 dipole pair, as well
as between the P1-P2 pair. This means that this array is most sensitive to resistivity changes
below the electrodes in each dipole pair. As the "n" factor is increased, the high sensitivity
values become increasingly more concentrated beneath the C1-C2 and P1-P2 dipoles, while
the sensitivity values beneath the center of the array between the C1-P1 electrodes decreases.
For "n" values of greater than 2, the sensitivity values at the pseudosection plotting point
become negligible. The sensitivity contour pattern becomes almost vertical for "n" values
greater than 2. Thus the dipole-dipole array is very sensitive to horizontal changes in resistivity,
but relatively insensitive to vertical changes in the resistivity. Thus it is good in mapping
vertical structures, such as dykes and cavities, but relatively poor in mapping horizontal
structures such as sills or sedimentary layers. The median depth of investigation of this array
depends on both the “a” spacing and the “n” factor (Table 2). In general, this array has a
shallower depth of investigation compared to the Wenner array, for example at n=1 the depth
of investigation is 0.416a compared to 0.512a for the Wenner Alpha array.
Figure 24. 2-D sensitivity sections for the dipole-dipole array. The sections with (a) n=1, (b)
n=2, (c) n=4 and (d) n=6.
Due to the almost vertical pattern of the sensitivity contours, the depth of investigation
(which a 1-D horizontal average of the sensitivity values) is not particularly meaningful for the
dipole-dipole array for "n" values greater than 2. From experience with synthetic modeling and
field data, the median depth of investigation might underestimate the depth of structures sensed
by this array by about 20% to 30% for the large “n” factors. For 2-D surveys, this array has
better horizontal data coverage than the Wenner (Figure 20). This can be an important
advantage when the number of nodes available with the multi-electrode system is small.
One possible disadvantage of this array is the very low signal strength for large values
of the “n” factor. The voltage is inversely proportional to the cube of the “n” factor. For the
same current, the voltage measured by the resistivity meter drops by about 56 times when “n”
is increased from 1 to 6 (Table 2). One method to overcome this problem is to increase the “a”
spacing between the C1-C2 (and P1-P2) dipole pair to reduce the drop in the potential when
the overall length of the array is increased to increase the depth of investigation. Figure 25
shows two different arrangements for the dipole-dipole array with the same array length but
with different “a” and “n” factors. The signal strength of the array with the smaller “n” factor
is about 28 times stronger than the one with the larger “n” factor.
To use this array effectively, the resistivity meter should have comparatively high
sensitivity and very good noise rejection circuitry, and there should be good contact between
the electrodes and the ground. With the proper field equipment and survey techniques, this
array has been successfully used in many areas to detect structures such as cavities where the
good horizontal resolution of this array is a major advantage.
The plotting location of the corresponding data point (based on the median depth of
investigation) used in drawing the apparent resistivity pseudosection is also shown in Figure
24. Note that the pseudosection plotting point falls in an area with very low sensitivity values
for “n” values of 4 and above. For the dipole-dipole array, the regions with the high sensitivity
values are concentrated below the C1-C2 electrodes pair and below the P1-P2 electrodes pair.
In effect, the dipole-dipole array gives minimal information about the resistivity of the region
surrounding the plotting point, and the distribution of the data points in the pseudosection plot
does not reflect the subsurface area mapped by the apparent resistivity measurements. Note
that if the data point is plotted at the point of intersection of the two 45 angle lines drawn from
the center of the two dipoles, it would be located at a depth of 0.7 units in Figure 24d (compared
with 0.19 units given by the median depth of investigation method) where the sensitivity values
are almost zero!
Loke and Barker (1996a) used an inversion model where the arrangement of the model
blocks directly follows the arrangement of the pseudosection plotting points. This approach
gives satisfactory results for the Wenner and Wenner-Schlumberger arrays where the
pseudosection point falls in an area with high sensitivity values (Figure 23a and Figure 26).
However, it is not suitable for arrays such as the dipole-dipole and pole-dipole where the
pseudosection point falls in an area with very low sensitivity values. The RES2DINV program
uses a more sophisticated method to generate the inversion model where the arrangement the
model blocks is not tightly bound to the pseudosection.
A final minor note. In most textbooks, the electrodes for this array are arranged in a C1-
C2-P1-P2 order that will in fact give a negative geometric factor. The arrangement assumed in
these notes is the C2-C1-P1-P2 arrangement.
Figure 25. Two possible different arrangements for a dipole-dipole array measurement. The
two arrangements have the same array length but different “a” and “n” factors resulting in
very different signal strengths.
Figure 26. 2-D sensitivity sections for the Wenner-Schlumberger array. The sensitivity
sections with (a) n=1, (b) n=2, (c) n=4 and (d) n=6.
Figure 27. A comparison of the (i) electrode arrangement and (ii) pseudosection data pattern
for the Wenner and Wenner-Schlumberger arrays.
This array has the widest horizontal coverage and the deepest depth of investigation.
However, it has the poorest resolution, which is reflected by the comparatively large spacing
between the contours in the sensitivity function plot (Figure 28).
𝐼𝜌
𝑉𝐶2 = 2𝜋𝑘(𝑘+1)𝑎.
The ratio of the potentials due to the C2 and the C1 electrodes is then given by
𝑉𝐶2 𝑛(𝑛+1) 𝑛 2
= ≈ (𝑘 )
𝑉𝐶1 𝑘(𝑘+1)
if n is reasonably large (example 10). Thus the effect of the C2 electrodes decreases with square
of its distance from the survey line. The pole-dipole array measurements are less affected by
the C2 remote electrode compared to the pole-pole array. If the distance of the C2 electrode is
more than 5 times the largest C1-P1 distance used, the error caused by neglecting the effect of
the C2 electrode is less than 5% (the exact error also depends on the location of the P2 electrode
for the particular measurement and the subsurface resistivity distribution).
Due to its good horizontal coverage, this is an attractive array for multi-electrode
resistivity meter systems with a relatively small number of nodes. The signal strength is lower
compared with the Wenner and Wenner-Schlumberger arrays but higher than the dipole-dipole
array. For I.P. surveys, the higher signal strength (compared with the dipole-dipole array)
combined with the lower EM coupling (compared with the Wenner and Wenner-Schlumberger
arrays) due to the separation of the circuitry of the current and potential electrodes makes this
array an attractive alternative.
The signal strength for the pole-dipole array decreases with the square of the “n” factor.
While this effect is not as severe as the dipole-dipole array, it is usually not advisable to use
“n” values of greater than 8 to 10. Beyond this, the “a” spacing between the P1-P2 dipole pair
should be increased to obtain a stronger signal strength. The ratio of the depth of investigation
to the array length does not increase significantly for ‘n’ values beyond 4 or 5 (Table 2). Thus
to increase the depth of investigation, it is better to increase the ‘a’ spacing rather than to keep
increasing the ‘n’ value when it reaches 5. There is another interesting effect when the ‘n’ factor
is increased that is frequently not appreciated, and leads to an interesting pitfall in field surveys.
This is discussed in section 4.10.
Figure 30. The pole-dipole array 2-D sensitivity sections. The sensitivity sections with (a)
n=1, (b) n=2, (c) n=4 and (d) n=6.
The value of xp where this integral reaches half the total value can then be used.
Figure 31. Sensitivity sections for the gradient array. The current electrodes are fixed at x=0
and 1.0 meters, and the distance between the potential electrodes is 0.1 m. The position of the
P1-P2 potential electrodes at (a) 0.45 and 0.55 m., (b) 0.55 and 0.65 m., (c) 0.65 and 0.75 m,
(d) 0.75 and 0.85 m. and (e) 0.80 and 0.90 m.
Figure 32. (a) Example of multiple gradient array data set and inversion model. (b) Profile
plot using exact pseudodepths. (c) Profile plot using approximate pseudodepths.
Figure 33. The apparent resistivity pseudosection for the dipole–dipole array using
overlapping data levels over a rectangular prism. Values of 1 to 3 meters are used for the dipole
length ‘a’, and the dipole separation factor ‘n’ varies from 1 to 5. Compare this with Figure 20c
for the same model but with ‘a’ fixed at 1 meter, and “n” varying from 1 to 10. In practice, a
‘n’ value greater than 8 would result in very noisy apparent resistivity values.
1 i 8
FI F2 Di x, z
8 i 1
(2.7)
Note that the absolute value of the sensitivity value is used, since negative sensitivity values
also give information about the subsurface. This is also to avoid a situation where the negative
sensitivity value for one measurement cancels out the positive value for another measurement.
Figure 34a shows a dipole-dipole array based measurement sequence. The two current
electrodes are set at one end of the line and the potential measurements are made using
successive pairs of the potential nodes. Note the area with the highest cumulative sensitivity is
around the first 4 electrodes. The area of deepest penetration is between the second and third
electrodes, i.e. between the current dipole and the first potential electrode.
The second configuration places the current electrodes at the ends of the line (Figure
34b). Almost all the measurements are made using a fixed separation of one unit electrode
spacing between the potential electrodes. The potential dipole is moved from one end of the
survey line to the other end, i.e. a gradient array type of arrangement. One advantage of this
configuration over the dipole-dipole arrangement is the larger voltage (i.e. lower noise level)
measured by the potential electrodes. The cumulative sensitivity section shows a more uniform
pattern compared to that produced by the dipole-dipole configuration. The depth of
investigation is slightly deeper below current electrodes, and slightly less below the center of
the line.
The third configuration also places the current electrodes at the ends of the lines, but
keeps the center of the potential dipoles at or near the center of the line. The first 4
measurements use a symmetrical Wenner-Schlumberger arrangement and increase the
separation between the potential electrodes until they reach next to the current electrodes. The
second set of 4 measurements uses separations of 2 and 3 times the unit spacing for the potential
dipole pair but with a mid-point slightly to one side of the center of the line. The cumulative
sensitivity pattern has an even more uniform pattern below the line compared with the gradient
array configuration.
To get an idea of the relative depths of investigation for the 3 configurations, we use
the light blue contour (with a value of 8 units) as a guide. The expanding Wenner-Schlumberger
configuration has the deepest depth with high sensitivity values while the dipole-dipole
configuration appears to have the shallowest. Note the current electrodes for dipole-dipole
array is limited to the left end of the line, and for this reason the region with the deepest depth
of investigation is skewed towards the left.
One problem with using the cumulative sensitivity sections is that it is a rather crude
measure of the resolution that the data set has. It probably overestimates the depth of
investigation for the gradient and Wenner-Schlumberger configurations (Figure 34b, c). This
is because it does not take into account the orthogonality of data, i.e. the information from two
measurements for the same region might not be independent. A better method to evaluate the
information content from a group of array configurations is the model resolution (section
4.9.6).
Figure 34. Cumulative sensitivity sections for different measurement configurations using (a)
a dipole-dipole sequence, (b) a moving gradient array and (c) an expanding Wenner-
Schlumberger array.
Figure 35. The output from the RES2DMOD software for the SINGLE_BLOCK.MOD 2-D
model file. The individual cells in the model are shown in the lower figure, while the upper
figure shows the pseudosection for the Wenner Beta (dipole-dipole with n=1) array.
cursor to one of the color boxes in the legend above the model, and click the resistivity value
you want. Press the F1 key to get information about the keys used by the program to edit the
model. Note that clicking the cells with the mouse buttons will only change the resistivity of
the cells displayed on the screen, but will not change the resistivity of the buffer cells towards
the left, right and bottom edges of the mesh. To change the resistivity of the buffer cells, you
need to use the “[“, “]” and “D” keys.
Next select the “Model Computation” option to calculate the potential values for this
model. The calculations will probably only take a few seconds, after which you should go back
to “Edit/Display Model” option. In this option, select the “Edit model” sub-option again to see
the apparent resistivity pseudosection for this model. The program will ask you to select the
type of contour intervals you wish to use. For most resistivity pseudosections choose the
‘Logarithmic contour intervals’, while for I.P. pseudosections choose the linear contour
intervals.
To change the type of array, use the “Change Settings” sub-option. Select another array,
such as the pole-pole or dipole-dipole, and see what happens to the shape of the contours in the
pseudosection.
the thickness and/or width of the prism and see what happens to the
anomaly in the pseudosection. Also try changing the resistivity of
the prism to 500 m.
BLOCK_TWO2.MOD (1). The default array type is the Wenner Alpha. Run the ‘Model
In this case, we have two computation’ option to get things started.
prisms at a mean depth of 1 (2). Next, take a look at the pseudosection. Is it possible to separate
meter, and 2 meters apart. The the highs due to each prism?
two prisms are identical and at (3). Try the whole range of arrays, i.e. the 3 different versions of
the same depth. We will try the Wenner, pole-dipole, dipole-dipole, Schlumberger etc. Which
different arrays to which ones arrays are more likely to be able to resolve the two prisms?
are more likely to be able to
resolve the prisms in the
pseudosection.
BLOCKS_UP.MOD (1). The default array type is the Wenner Alpha. Run the ‘Model
Now we have two prisms, one computation’ option to get things started.
on top of another. Is it possible (2). Next, take a look at the pseudosection. Is it possible to separate
to tell them apart? the highs due to each prism?
(3). Try the whole range of arrays, i.e. the 3 different versions of
the Wenner, pole-dipole, dipole-dipole, Schlumberger etc. Which
arrays are more likely to be able to resolve the two prisms?
THICK_DIKE.MOD (1). The default array type is the Wenner Beta. Calculate the
This model has a wide low potentials and take a look at the pseudosection. Is it possible to infer
resistivity dike in high the existence of the dike from the pseudosection?
resistivity bedrock. (2). Change to the Wenner Alpha, and see what happens.
(3). Try other arrays such as the Wenner Gamma, Wenner-
Schlumberger and pole-dipole. Which arrays show a significant
low resistivity anomaly?
THIN_DIKE.MOD (1). The default array type is the Wenner Beta. Calculate the
This model has a thin dike in potentials and take a look at the pseudosection. Is it possible to infer
high resistivity bedrock. This is the existence of the dike for the pseudosection?
a more difficult target than the (2). Change to the Wenner Alpha, and see what happens.
thick dike. (3). Try other arrays such as the Wenner Gamma, Wenner-
Schlumberger and pole-dipole. Which arrays show a significant
low resistivity anomaly?
MODEILIP.MOD (1). The default array type is the dipole-dipole, which is probably
An I.P. model just to round the most commonly used array for I.P. surveys. Calculate the
things up. potential values, and take a look at the pseudosections.
(2). Another array that is sometimes used for I.P. surveys is the
pole-dipole that has the advantage of a stronger signal strength.
Check the pseudosections obtained with this array.
4.1 Introduction
After the field survey, the resistance measurements are usually reduced to apparent
resistivity values. Practically all commercial multi-electrode systems come with the computer
software to carry out this conversion. In this section, we will look at the steps involved in
converting the apparent resistivity values into a resistivity model section that can be used for
geological interpretation. I will assume that the data is already in the RES2DINV format. The
conversion program is provided together with many commercial systems. So far, the ones that
have the conversion program include Abem, AGI, Campus, Geofysika, Geometrics, Iris, OYO,
Pasi and Scintrex. If your equipment manufacturer is not on the list, please contact them about
the conversion software. The data format used by the RES2DINV program is described in detail
in the RES2DINV.PDF manual provided with the program. Please refer to the manual for the
details.
In the next section, we will look at two methods to handle bad data points. Such bad
data points should be removed before a final interpretation is made.
Due to the large variety of data sets collected over various geological environments, no
single inversion method will give the optimum results in all cases. Thus the RES2DINV
program has a number of settings that can be changed by the user to obtain results that are
closer to the known geology. The various options are also discussed in the following sections.
can be used for practically any array and any distribution of the data points. The main
disadvantage of the method is the larger amount of computer time needed. In this method, a
preliminary inversion of the data set is first carried with all the data points. After carrying out
the trial inversion, switch to the 'Display' window in RES2DINV, and read in the INV file
containing the inversion results. After that, select the 'RMS error statistics' option that displays
the distribution of the percentage difference between the logarithms of the measured and
calculated apparent resistivity values. The error distribution is shown in the form of a bar chart,
such as in Figure 38. Normally, the highest bar is the one with the smallest errors, and the
heights of the bars should decrease gradually with increasing error values. The bad data points,
caused by problems such as poor ground contact at a small number of electrodes, should have
significantly higher errors than the "good" data points.
Figure 36. An example of a field data set with a few bad data points. The most obvious bad
datum points are located below the 300 meters and 470 meters marks. The apparent resistivity
data in (a) pseudosection form and in (b) profile form.
Figure 37. Selecting the menu option to remove bad data points manually.
Figure 38 shows the error distribution bar chart for the data set shown in Figure 36 that
has a few bad data points. In the bar chart, almost all the data points have errors of 20 percent
or less. The bad data points show up data points with errors of 60 percent and above, which
can be easily removed from the data set by moving the green cursor line to the left of the 60%
error bar. In this way the 5 bad data points are removed from this data set. For some data sets,
the error distribution might show a more complicated pattern. As a general rule, data points
with misfits of 100 percent and above can usually be removed.
Figure 38. Error distribution bar chart from a trial inversion of the Grundfor Line 1 data set
with five bad data points.
given by
J T J W T W ΔqJ T g W T Wq , (4.2)
so that the model resistivity values, q, changes in a smooth manner. The next option ‘Use
combined inversion method’ attempts to combined the smoothness-constrained method as
given in (4.1) with the Marquardt-Levemberg as given in (1.20). However, the result obtained
by this combination has not been very impressive and will not be examined.
The ‘Select robust inversion’ option has proved to be much more useful. It modifies the
formulation in (4.2) so that different elements of the model parameter change and data misfit
vectors have the same magnitudes. It is given by
J T R d J W T R m W Δqk J T R d g W T R m Wq k . (4.3)
The details are described in section 1.4. This method is also known as an l 1-norm or
robust or blocky inversion method (Loke et al. 2003), whereas the conventional smoothness-
constrained least-squares method as given in equation (4.2) is an l 2-norm inversion method.
The l2-norm inversion method gives optimal results where the subsurface geology exhibits a
smooth variation, such as the diffusion boundary of a chemical plume. However, in cases where
the subsurface consists of bodies that internally homogeneous with sharp boundaries (such as
an igneous dyke), this method tends to smear out the boundaries. The l 1-norm or blocky
optimization method tends to produce models that are piecewise constant (Ellis and Oldenburg,
1994a). This might be more consistent with the known geology in some situations.
Figure 40 shows the inversion results for data from a synthetic model with sharp
boundaries. In this case, the robust inversion method gives significantly better results since the
true model consists of three homogenous regions with sharp boundaries. Many synthetic test
models are of a similar nature with sharp boundaries, so not surprisingly, results obtained with
the l2-norm smooth inversion method are not optimal for such data sets (Olayinka and
Yaramanci, 2000).
Resistivity values have a logarithmic range, possibly ranging from less than 1 to over
1000 m in a single data set. By using the logarithm of the resistivity values as the parameters
in the inversion process, the numerical range of the parameters can be reduced to a linear range.
However, in some situations, it is not possible to make use of the logarithm if there are negative
or zero values. This does not usually occur for normal surface surveys with the standard arrays,
but could occur in borehole surveys or if non-standard arrays are used. The program also allows
the use to use the apparent resistivity value directly as the inversion parameter in the “Choose
Figure 40. Example of inversion results using the l2-norm smooth inversion and l1-norm blocky
inversion model constrains. (a) Apparent resistivity pseudosection (Wenner array) for a
synthetic test model with a faulted block (100 m) in the bottom-left side and a small
rectangular block (2 m) on the right side with a surrounding medium of 10 m. The
inversion models produced by (b) the conventional least-squares smoothness-constrained or l2-
norm inversion method and (c) the robust or l1-norm inversion method.
The setting ‘Type of method to solve least-squares equation’ provides two numerical methods
to solve the least-squares equation (as in equations 4.1 to 4.3). The default method, particularly
when the number of model parameters n is small, is the standard (or complete) Gauss-Newton
method where a direct method (Golub and van Loan, 1989) is used to solve the equation. This
method produces an exact solution but as the time taken is proportional to n3, it could take a
very long time for large data set with over 5000 model cells. An alternative method, the
incomplete Gauss-Newton method, can be used for such cases. In the incomplete Gauss-
Newton method, an approximate solution of the least method is determined by using an
iterative method (Golub and van Loan, 1989). The final solution obtained with this method has
a difference of about 1 to 2 percent compared to the complete Gauss-Newton method solution.
By sacrificing a small amount of accuracy in the solution to the least-squares equation, the
computer time required could be reduced by a factor of 5 to 10 times for very large models.
of data points. The thickness of each deeper layer is increased to reflect the decreasing
resolution of the resistivity method with increasing depth. In general, this produces a model
where the thickness of the layers increases with depth, and with thicker cells at the sides and
in the deeper layers (Figure 41a). For most cases, this gives an acceptable compromise. The
distribution of the data points in the pseudosection is used as a rough guide in allocating the
model cells, but the model section does not rigidly follow the pseudosection. However, the user
can change the width and thickness of the cells using a variety of options.
Figure 41. Different methods to subdivide the subsurface into rectangular prisms in a 2-D
model. Models obtained with (a) the default algorithm, (b) by allowing the number of model
cells to exceed the number of data points, (c) a model which extends to the edges of the
survey line and (d) using the sensitivity values for a homogeneous earth model.
After reading a data file, clicking the ‘Display model blocks’ option will show the
distribution of the model cells currently used by the program, such as in Figure 41. Clicking
the ‘Change thickness of layers’ option will bring up the following dialog box.
Figure 42. The options to change the thickness of the model layers.
When the program reads in a data file, it will normally set the first layer thickness using
the minimum pseudodepth of the data points. For surface surveys, since the resolution
decreases with depth, the thickness of the layers is normally increased by between 5 to 15 %
with each deeper layer. The program normally uses a model where the depth to the deepest
layer does not exceed the maximum pseudodepth in the data set. To use a model that spans a
deeper depth range, you can change the factor to increase model depth range, for example from
1.0 to 1.30 to increase the model depth range by 30%. In this dialog box, you can also allow
the program to use a model where the number of model cells exceeds the number of data points.
This is useful to avoid having a model with very wide cells near the bottom for data sets with
very sparse data sets. Figure 41b shows such an example with this option enabled.
The ‘Modify depth to layers’ option allows you to set the depth to each layer
individually. This is useful if you want a layer boundary to coincide exactly with a known
geological boundary.
In ‘Use extended model’ option, model cells of uniform thickness right up to the left
and right edges of the survey line are used (Figure 41c). This is probably an extreme case. As
the number of model parameters increase, the computer time needed to carry out the inversion
also increases. This can be an important consideration for very large data sets with several
hundred electrodes. The ‘Make sure model blocks have same widths’ option is probably more
useful. It uses a base model such as in Figure 41b, but avoids thicker cells at the sides. It ensures
that the cells at the sides to have the equal widths.
The ‘Reduce effect of side blocks’ option affects the calculation of the Jacobian matrix
values for the model cells located at the sides and bottom of the model section. Normally, for
a cell located at the side, the contributions by all the mesh elements associated with the model
cell are added up right to the edge of the mesh. This gives a greater weight to the side cell
compared to the interior cells. In some cases, particularly when the robust inversion option is
used, this can result in unusually a high or low resistivity value for the side cell. This option
leaves out the contribution of the mesh elements outside the limits of the survey line to the
Jacobian matrix values for the side cells.
The last two options, ’Change width of blocks’ and ‘Type of cross-borehole model’,
are only used for certain special cases, and will not be described here. Please refer to the
RES2DINV.PDF manual for the details.
c). Model sensitivity options
This set of options relate to the model sensitivity (Jacobian matrix) values. The ‘Display
model blocks sensitivity’ option will display the sum of the absolute values of the sensitivity
values associated with the model cell. Figure 41d shows an example. Note the blocks at the
sides and bottom has greater sensitivity values due to the larger sizes. To avoid this effect, the
‘Display subsurface sensitivity’ option divides the subsurface into cells of equal size and shows
the sensitivity values. This is useful to get an idea of the regions of the subsurface “scanned”
by the survey configuration used.
All the techniques used to subdivide the subsurface described in the earlier section are
based on heuristic algorithms. Figure 41d shows the block distribution generated by a more
quantitative approach based the sensitivity values of the model cells. This method is selected
by the ‘Generate model block’ option. This technique takes into account the information
contained in the data set concerning the resistivity of the subsurface for a homogeneous earth
model. It tries to ensure that the data sensitivity of any cell does not become too small (in which
case the data set does not have much information about the resistivity of the cell).
The next set of inversion options are grouped below the ‘Change Settings’ choice on
the main menu bar. Clicking this will bring up the list of options shown below.
Figure 43. The options under the ‘Change Settings’ menu selection.
Since the model resolution decreases with depth, the program increases the damping
factor value by about 5% for each deeper layer. You can change this factor using the “Change
damping factor with depth’ option. Instead of automatically decreasing the damping factor by
half after each iteration, the ‘Optimize damping factor’ option allows the program to look for
an optimum damping factor. The time taken per iteration will be more, but fewer iterations
might be needed for the program to converge.
The ‘Limit range of model resistivity values’ option is intended for cases where the
default settings (without limits) produces a model with resistivity values that are too high or
too low. Selecting this option will bring up the dialog box shown in Figure 44. In this example,
the upper limit for is 20 times the average model resistivity value for the previous iteration
while the lower limit is 0.05 times (i.e. 1/20 times). The program uses “soft” limits that allow
the actual resistivity model values to exceed the limits to a certain degree. However, this option
will avoid extremely small or large model resistivity values that are physically unrealistic.
The ‘Vertical/Horizontal flatness filter ratio’ option allows the user to fine-tune the
smoothness-constrain to emphasize vertical or horizontal structures in the inversion model.
Some geological bodies have a predominantly horizontal orientation (for example sedimentary
layers and sills) while others might have a vertical orientation (such as dykes and faults). This
information can be incorporated into the inversion process by setting the relative weights given
to the horizontal and vertical flatness filters. If for example the structure has a predominantly
vertical orientation, such as a dyke, the vertical flatness filter is given a greater weight than the
horizontal filter.
Figure 44. The dialog box to limit the model resistivity values.
equation 2.17 in section 2.5.12) is a measure of the amount of information about the resistivity
of a model block cell in the measured data set. The higher the sensitivity value, the more reliable
is the model resistivity value. In general, the cells near the surface usually have higher
sensitivity values because the sensitivity function has very large values near the electrodes. The
cells at the sides and bottom can also have high sensitivity values due to the much larger size
of these cells that are extended to the edges of the finite-difference or finite-element mesh (the
program has an option to reduce this effect which might produce artifacts at the edges of the
model). If you had carried out an inversion of the data set before calling this option, the program
will make use of the Jacobian matrix of the last iteration. Otherwise, it will calculate the
Jacobian matrix for a homogenous earth model. Figure 41d shows an example of a plot of the
sensitivity section for a model.
Figure 45b shows the model section obtained from the inversion of a data set for a
survey to map leakage of pollutants from a landfill site (Niederleithinger, 1994). The model
sensitivity section in Figure 45c shows high sensitivity values near the surface with decreasing
values with depth. This is to be expected as the near surface materials have a larger influence
on the measured apparent resistivity values.
Figure 45. Landfill survey example (Wenner array). (a) Apparent resistivity pseudosection ,
(b) model section, (c) model sensitivity section (d) model uncertainty section, (e) minimum
and maximum resistivity sections.
In order to assess the accuracy of the inversion model, an estimate of the reliability of the model
values is required. One possible approach is by using the model covariance matrix (Menke,
1984). This is commonly used for models that consist of a small number of parameters (such
as the 1-D model in Figure 7). Figure 45d shows the model uncertainty values obtained from
the covariance matrix method as described by Alumbaugh and Newman (2000) where the
smoothness constraint is included in the model uncertainty estimate. In this way, the model
uncertainty values are less sensitive to size of the model cells. However, the uncertainty values
are only meaningful if the subsurface resistivity varies in a smooth manner, as assumed by the
smoothness constraint. If the subsurface resistivity does not vary in a smooth manner, this
method is likely to underestimate the actual uncertainty. Figure 45e shows the maximum and
minimum resistivity values of each cell at the limits of the model uncertainty range. Features
that are common to both model sections can be considered more reliable.
J
J FR Δqk J T R d g FR (q k q o ) ,
T
(4.4)
with FR s x C x R m C x z C z R m C z .
T T
Figure 46. Landfill survey depth of investigation determination. (a) Model section with
extended depths and (b) the normalized DOI index section.
Figure 47. Beach survey example in Denmark. The figure shows the model section with
extended depths and the normalized DOI index section.
The DOI method is useful in marking the regions where the model values are well
constrained by the data set, and thus greater confidence can be placed on the model resistivity
values at such regions. The DOI method may be considered an empirical method to determine
the regions where we can reasonably resolve the subsurface. Another method based on the
model resolution values is discussed in section 4.9.
Figure 48. Different methods to incorporate topography into a 2-D inversion model. (a)
Schematic diagram of a typical 2-D inversion model with no topography. A finite-element
mesh with four nodes in the horizontal direction between adjacent electrodes is normally used.
The near surface layers are also subdivided vertically by several mesh lines. Models with a
distorted grid to match the actual topography where (b) the subsurface nodes are shifted
vertically by the same amount as the surface nodes, (c) the shift in the subsurface nodes are
gradually reduced with depth or (d) rapidly reduced with depth, and (e) the model obtained
with the inverse Schwartz-Christoffel transformation method.
GLADOE2.DAT – A data set (1). If you have the time, run the same tests as you did
taken to check for leakage from earlier for the RATHCRO.DAT data set.
a dam (Dahlin pers. comm.). This is another example where the option to reduce the
The survey area has effect of the side cells makes a significant difference when
topography. the robust inversion option is used.
The first item is the number of regions where the resistivity is to be specified. In the example
above, 2 regions are specified. If a value of 0 is given (default value), then there is no region
where the resistivity is specified by the user. Next, the shape of the region is given, R for
rectangular or T for triangular. If a rectangular region is specified, then the X and Z coordinates
of the top-left and bottom-right corners of the rectangle are given, as shown in the Figure 49.
If a triangular region is chosen, the X and Z coordinates of the 3 vertices of the triangle must
be given in an anti-clockwise order. After the coordinates of the region to be fixed are given,
the next data item is the resistivity of the region. After that, the damping factor weight for the
resistivity of the region is needed. This parameter allows you control the degree in which the
inversion subroutine can change the resistivity of the region. There is usually some degree of
uncertainty in resistivity of the region. Thus, it is advisable that the program should be allowed
(within limits) to change the resistivity of the region. If a damping factor weight of 1.0 is used,
the resistivity of the region is allowed to change to the same extent as other regions of the
subsurface model. The larger the damping factor weight is used, the smaller is the change that
is allowed in the resistivity of the "fixed" region. Normally, a value of about 1.5 to 2.5 is used.
If a relatively large value is used, for example 10.0, the change in the resistivity of the region
would be very small during the inversion process. Such a large value should only be used if the
resistivity and shape of the region is accurately known. Figure 50 shows the allocation of the
cells in the subsurface model together with the fixed regions for the MODELFIX.DAT data
set.
Seismic refraction and seismic reflection surveys are commonly used in engineering
surveys. Both methods can give accurate and detailed profiles of the subsurface interfaces. In
some cases, a distinct and sharp transition between two layers can be mapped by the seismic
survey. This information can be used to improve the results from the inversion of a 2-D
resistivity imaging survey along the same line. The subsurface in the inversion model can be
divided into two zones, one above and one below the interface calculated from the seismic
survey. The resistivity values are constrained to vary in a smooth manner within each zone, but
an abrupt transition across the zone boundary is allowed by removing any constrain between
the resistivity values below and above the zone boundary (Smith et al. 1999).
Figure 51 shows an example where the boundaries of a clay and a gravel layer were
known from a seismic refraction survey (Scott et al., 2000). The known boundaries were then
incorporated into the resistivity inversion model (Figure 51a) that allows the model resistivity
to change abruptly across the boundaries. Note the sharp contrast across the gravel layer and
the underlying low resistivity clay layer (marked by the bottom blue region in lower section in
Figure 51b).
Figure 49. Fixing the resistivity of rectangular and triangular regions of the inversion model.
Figure 50. The inversion model cells with fixed regions. The fixed regions are drawn in purple.
Note that the triangular region extends beyond the survey line.
Figure 51. Example of an inversion model with specified sharp boundaries. (a) The boundaries
in the Clifton survey (Scott et al., 2000) data set is shown by the blue lines. (b) The measured
apparent resistivity pseudosection and the inversion model section.
Figure 52. The effect of cell size on the model misfit for near surface inhomogeneities. (a)
Model with a cell width of one unit electrode spacing. (b) A finer model with a cell width of
half the unit electrode spacing. The near surface inhomogeneities are represented by coloured
ovals.
Figure 53 shows a synthetic model used to illustrate the effect of the model cell misfit.
The main structure is a faulted block of 100 m and a rectangular prism of 1 m in a medium
of 10 m. A series of small near-surface high resistivity blocks with widths of 1.0, 0.75, 0.50
and 0.25 m. and resistivity of 300 m are placed above the faulted block. A similar series of
low resistivity blocks of 1.0 m are located to the left.The pole-dipole array has the P1-P2
dipole length (“a”) fixed at 1.0 m., but with “n” factor ranging from 1 to 16. Note the strong
anomalies produced by the near-surface inhomogeneities. Note also that the Wenner array is
much less affected by the near-surface inhomogeneities. The reason lies in the sensitivity
patterns of the two arrays (compare Figure 23 and Figure 30). For the pole-dipole array with
large “n” values, the region with the highest positive sensitivity values is concentrated below
the P1-P2 dipole pair.
Figure 54 shows the inversion results for the pole-dipole data set with different cell
widths. The model with a cell width of 1.0 m. (i.e. one unit electrode spacing) shows significant
distortions near the top of the faulted block as well as in the low resistivity rectangular block
(Figure 54b). Most of the distortions are removed in the model with a cell width of half the unit
electrode spacing (Figure 54c). This means the residual misfits with widths of up to one-quarter
the unit electrode spacing does not have a significant effect on the calculated apparent
resistivity values. The model with a cell width of one-quarter the unit electrode spacing (Figure
54d) does not show any major improvement over the half-cell width model although in theory
it should more accurately model the near surface inhomogeneities. In fact, there is a poorer
agreement with the true model in the lower part of the faulted block (Figure 54d). Experiments
with a number of data sets show that using a cell width of one-quarter the unit electrode spacing
can sometimes lead to oscillating model resistivity values. This is probably because the data
does not have sufficient information to accurately resolve such small cells.
Figure 53. Synthetic model (c) used to generate test apparent resistivity data for the pole-dipole
(a) and Wenner (b) arrays.
Figure 55 shows an example from a survey over an underground pipe using the Wenner-
Schlumberger array. There are large resistivity variations near the surface, probably due to
stones in the topmost layer of the soil. If the effect of the near surface variations are not
accurately accounted for by the inversion model, it can lead to distortions in the lower portions
of the model as the programs attempts to reduce the data misfit by distorting the lower part of
the model. The model with a cell width of half the unit electrode spacing in Figure 55c shows
a slightly better fit with the measured data (i.e. lower RMS error) and a more circular shape for
the low resistivity anomaly below the 12 meters mark.
In conclusion, for most cases, using a cell width of half the unit electrode spacing seems
to give the optimum results. Using a cell width of one-third the unit spacing seems to be
beneficial only a certain cases with the pole-dipole and dipole-dipole arrays with very large ‘n’
values (see section 4.10). A cell width of one-quarter the unit spacing sometimes leads to
instability with oscillating model values, particularly in the first few layers. Thus the use of a
cell width of less than one-quarter the true unit electrode spacing is not advisable.
Note that using finer cells will lead to longer inversion times, so using a width of the
half the unit electrode spacing seems to provide the best trade-off. These examples present a
strong case for using a model with a cell width of half the unit electrode spacing as the default
choice in the inversion of most data sets. It avoids the problem caused by model cells boundary
misfits if the cell is too coarse, and the increase in computer time is tolerable.
One effect of using finer model cells for some data set is the emergence of 'ripples' in
the first couple of layers in the inversion model. Figure 56a shows the pseudosection from a
survey in the Blue Ridge region in north-east USA (Seaton and Burbey, 2000). The inversion
model using the blocky model inversion norm in Figure 56b show some rapid near-surface
alternating resistivity values between the 160 and 240 meter marks. These artifacts can be
significantly reduced by using a higher damping factor for the top layer as shown in Figure
56c.
As a final note, it appears that the effect of using narrower model cells is less dramatic
in 3-D inversion. This is probably because in the 2-D model each cell extends to infinity in the
y-direction, whereas in the 3-D model the same cell is divided into a large number of much
smaller cells. Thus the effect of a single near-surface cell in the 3-D model on the calculated
apparent resistivity values is much smaller than in the equivalent 2-D model.
Figure 54. The effect of cell size on the pole-dipole array inversion model. Pole-dipole array.
(a) The apparent resistivity pseudosection. The inversion models obtained using cells widths
of (b) one, (b) one-half and (c) one-quarter the unit electrode spacing.
Figure 55. Example of the use of narrower model cells with the Wenner-Schlumberger array.
(a) The apparent resistivity pseudosection for the PIPESCHL.DAT data set. The inversion
models using (b) cells with a width of 1.0 meter that is the same as the actual unit electrode,
and (c) using narrower cells with a width of 0.5 meter.
Figure 56. Example of reduction of near-surface 'ripples' in inversion model. (a) Apparent
resistivity pseudosection for the BLUERIDGE.DAT data set. (b) Normal inversion using the
robust inversion norm and model refinement. (c) Inversion using higher damping factor for
first layer to reduce the 'ripple' effect in the top layer.
This will round up the positions of the electrodes for all the data points to the nearest unit
electrode spacing specified in the data file. Once the data in the correct format, we next look at
ways to reduce the computer time needed to invert the data set.
The program has options to use 2 or 4 nodes between adjacent electrodes. For a long survey
line, the program automatically selects the 2 nodes option. Here we will examine the effect of
using 2 or 4 nodes on the resulting inversion model. Normally selecting the 4 nodes option
increases the accuracy of the finite element routine used. However, due to the nature of the data
from a mobile survey, the effect might be minimal. There should be at least 4 horizontal nodes
between two electrodes in an array used in a measurement. In mobile surveys, the distance
between adjacent electrode positions is frequently much smaller than the electrode spacing in
an array. The example in Figure 59 has a spacing of 1 meter between adjacent electrode
positions but has an electrode spacing of 3 meter in an array. Thus there are 6 nodes between
two electrodes in the same array although a 2 nodes spacing is used between adjacent electrode
positions. As shown in Table 8, using 2 nodes instead of 4 nodes reduces the inversion time by
almost half (from 2093 to 1074 seconds) for the test data set. The effect on the resulting
inversion model is minimal as shown by Figure 57b,c. The difference in the data misfit is only
about 0.1% and there are no obvious differences in the model sections.
further reduces the calculation time by about 14% (Table 8), while the inversion model (Figure
57d) is essentially identical to the one obtained with the full Jacobian matrix (Figure 57c).
Table 8. Inversion times for the long Redas survey data set using different settings.
Trial Setting used Time (s)
1 Standard settings with 4 nodes between adjacent electrodes. 2093
2 Standard settings with 2 nodes between adjacent electrodes. 1074
3 Use fast calculation of Jacobian matrix values 934
4 Use sparse inversion techniques 210
5 Use 2 meter model cells 181
6 Use 3 meter model cells 171
Figure 57. Inversion models for a long survey line using different settings to reduce the
calculation time. (a) The measured apparent resistivity pseudosection for the Redas underwater
survey for the first 2600 meters. Inversion models using (b) standard inversion settings with 4
nodes between adjacent electrode positions, (c) standard inversion settings with 2 nodes
between adjacent electrode positions, (d) with calculation of Jacobian matrix values for
selected model cells, (e) with sparse inversion technique, (f) with 2 meter model cell width and
(g) with 3 meter model cell width.
Figure 58. Part of finite-element mesh used to model a survey with submerged electrodes. The
resistivity of mesh cells in the water layer are fixed at the known water resistivity.
Figure 59. Part of a finite-element mesh for a long survey line. The electrodes used in the same
array are 3 meters apart, while the data unit electrode spacing is 1 meter. Using 2 nodes between
adjacent electrodes in the mesh will actually result in 6 nodes between electrodes in the same
array.
(4.8)
If the cells are perfectly resolved, the model resolution matrix will have the form below
where the diagonal elements are 1.0 and elsewhere it is 0, so the model resolution matrix can
be written as
This means the calculated value for each cell only depends on the true value. In the
case with imperfect resolution, we might have something like
The diagonal elements give the ‘degree’ of resolution, while the off diagonal elements
give the degree of ‘contamination’ or cross-correlation with the neighboring model cells. One
way to illustrate the resolution is to plot the values of the diagonal elements of the R matrix.
This shows the degree at which the calculated model value depends on the true value. Some
authors choose a value of about 0.05 (5%) as the ‘cutoff’ value.
electrode system with 30 independent nodes and a spacing of 1.0 meter between adjacent
electrodes. First we look at the resolution for the widely used Wenner array. Assume we make
all the possible measurements, with the “a” spacing ranging from 1.0 to 9.0 meters giving a
total of 135 data points. The model resolution section (Figure 60a) shows that the resolution is
greatest near the surface, decreases rapidly with depth, and is very small below a depth of about
2.0 meters. Next we look at the Wenner-Schlumberger array with the ‘a’ spacing ranging from
1.0 to 3.0 m. and the ‘n’ factor ranging from 1 to 8 giving a total of 292 data points. It performs
significantly better than the Wenner array (Figure 60b) with significant resolution values up to
about 3.0 m. The third array we will look at is the dipole-dipole array with the dipole length
‘a’ fixed at 1.0 m. and the ‘n’ factor ranging from 1 to 6, giving a total of 147 data points. Its
performance is surprisingly good and comparable to the Wenner-Schlumberger array with
significant model resolution values to about 3.5 meters, although it only has about half the
number of data points (Figure 60c). It is much better than the Wenner array although it only
has 12 more data points.
In an attempt to improve the resolution for the dipole-dipole array, we next attempt to
use overlapping data levels with it. The ‘a’ dipole length ranges from 1.0 to 3.0 m. while the
‘n’ factor ranges from 1 to 6. The region with significant resolution values increases to about
5 meters (Figure 60d), an improvement from the 3.5 meters limit with a single ‘a’ dipole length
(although it is mainly concentrated near the center). It only has about 50 data points more than
the Wenner-Schlumberger (342 compared to 292), but it is a significant improvement.
It is possible to further improve the resolution we can get for a limited number of
independent nodes and data points?
Figure 60. Model resolution sections for the (a) Wenner, (b) Wenner-Schlumberger and (c)
simple dipole-dipole array (d) dipole-dipole array with overlapping data levels.
arrangement (with interleaved current and potential electrodes), and those with very high
geometric factors (i.e. low potentials) of greater than the dipole-dipole with a=1 and n=6. This
reduces the total number of possible array configurations to 51283, which is called the
‘comprehensive’ data set (Stummer et al., 2004; Wilkinson et al., 2006b). If we make
measurements with all these possible array configurations, we should get the best possible
resolution. The model resolution section for this ‘comprehensive’ data set is shown in Figure
61a. The region with significant resolution values extends to about 8 meters, significantly more
than the dipole-dipole (5 meters) array.
Figure 61. (a) Model resolution section for comprehensive data set. (b) Model resolution
section for optimized data set generated by the ‘Compare R’ method.
Although taking all the 51283 possible readings will give the maximum possible resolution,
this is not practical in an actual field survey. We want to take a small subset of the possible
readings that gives almost the same resolution. The main steps taken by the array optimization
routine is as follows.
1. Start with a small set, the ‘base’ data set, such as the dipole-dipole readings with a=1
and n=1 to 6 (147 readings).
2. Calculate the increase in the model resolution as each new array is added to the ‘base’
set (there are several alternatives for this step with differences in speed and accuracy).
Add the arrays that give the highest increase to the model resolution (and are
sufficiently independent of each other) to the base set. After each iteration, increase the
base set by a set percentage (usually about 1% to 5%).
3. Use the new set (base plus new arrays) as the base set, and repeat the procedure. Stop
when the desired maximum number of arrays (eg. 5000) is reached.
Different Methods have been proposed to carry out step (2) in selecting the arrays that
will give the maximum increase in the model resolution (Wilkinson et al., 2006b). Here, only
the results of one method, the ‘Compare R’ method, will be shown. It gives the most accurate
results in terms of optimized arrays, but it also takes the longest computer time to generate the
arrays (Wilkinson et al., 2006b). However the time needed has been reduced by more than
three orders of magnitude using optimized computer code and parallel programming techniques
(Loke et al., 2010a, 2010b, 2015a). The model resolution section for optimized data set (Figure
61b) is almost identical with that for the ‘comprehensive’ data set (Figure 61a). The optimized
data set has only 4462 array configurations, which is about 9% of the ‘comprehensive’ data set.
To illustrate the performance of the different array configurations in detecting
subsurface structures, a test model with four rectangular blocks at different depths (Wilkinson
et al., 2006b) is used (Figure 62). Synthetic data sets were generated for the Wenner, Wenner-
Schlumberger, dipole-dipole and optimized arrays. The data sets were then inverted to recover
back the subsurface resistivity.
Figure 63a shows the resulting inversion model for the Wenner array data set. The two
topmost blocks are well resolved while the third block is poorly resolved. The deepest block is
completely unresolved. The Wenner-Schlumberger array data set (Figure 63b) performs slight
better in resolving the third block but still cannot resolve the deepest block. The simple dipole-
dipole array data set with a single dipole length (a=1.0m.) performs significantly better than
even the Wenner-Schlumberger array data set where the third block is well resolved (Figure
63a), even though it has fewer data points (147 compared with 292). However, it still cannot
resolve the deepest block. The dipole-dipole array data set with multiple dipole lengths (all
possible data points with a=1.0 to 3.0 meters and n=1 to 6, with the restriction that the geometric
factor cannot exceed that for a=1.0 meter and n=6) shows higher resistivity values at the
location of the deepest block but it is not able to resolve the block (Figure 63d).
Figure 62. Test inversion model for the different arrays. (a) Apparent resistivity pseudosection
for the simple dipole-dipole array with a dipole length ‘a’ of 1.0 meter with the dipole
separation factor ‘n’ of 1 to 6. (b) Apparent resistivity pseudosection for the dipole-dipole array
with overlapping data levels. (c) The synthetic model with 4 rectangular blocks of 100 .m
embedded in a medium of 10 .m. as used by Wilkinson et al. (2006b).
The deepest block is fairly well resolved in the inversion model for the optimized data
set with 4462 data points (Figure 64a) where the shape and dimensions are close to the true
structure. It might be argued that the better performance of the optimized data set compared to
the dipole-dipole array with multiple dipole lengths is due to the much larger number of data
points (4462 compared to 342). Next we use a truncated subset of the optimized data set with
only 413 data points that is comparable to the dipole-dipole array data set. The inversion model
with the truncated optimized data set (Figure 64b) still shows significant better resolution for
the deepest block compared to the dipole-dipole array model (Figure 63d). While the shape of
the deepest block is significantly less sharp compared to the large optimized data set, it still is
better defined compared to the dipole-dipole array model. The third deepest block is also much
better resolved by both optimized data sets with resistivity values of over 50 Ohm.m at the
center compared to less than 30 Ohm.m for the dipole-dipole array model.
This simple example shows that there can be significant improvements in the depth
resolution by using an optimized data set. An old rule of thumb in resistivity sounding is that
the maximum depth of investigation is about one-sixth (17%) the line length, i.e. slightly less
than 6 meters near the center of the line in the above examples. However, using an optimized
data set has pushed this limit close to 8 meters or slightly more than one-fifth (about 28%) the
line length.
The array optimization problem for 2-D surveys has been quite well studied (Wilkinson
et al., 2012). Loke et al. (2014) describes the use of array optimization for cross-borehole
surveys. With recent advances in computer software and hardware, it can be used for long 2-D
survey lines (Loke et al., 2015a) and even 3-D surveys (Loke et al., 2014c). It is interesting to
note since the initial work by Wilkinson et al. (2006b), the calculation time required by the
‘Compare R’ method to generate the optimized arrays for a survey line with 30 electrodes has
been reduced from about 6 hours to 3 seconds, an improvement of about 7,000 times (Loke et
al., 2015a) over 9 years!
Figure 63. Inversion results with the (a) Wenner array, (b) Wenner-Schlumberger array, (c)
simple dipole-dipole array and (d) dipole-dipole array with overlapping data levels.
Figure 64. Inversions models with the (a) optimized data set with 4462 data points and (b) a
truncated optimized data set with 413 data points.
4.9.5 DOI versus model resolution?
In section 4.4 we have looked at the Depth of Investigation (DOI) method to estimate
the practical depth of investigation of field data sets. How does the model resolution section
compare with the DOI? The model sections and DOI section for the Landfill survey data set
were shown in Figure 45 and Figure 46. Figure 65 shows a comparison between the model
resolution and DOI sections for this data set. Using a value of 0.05 as the cutoff point for the
model resolution values might be arbitrary as it does not take into account the model
discretization used. If a finer model discretization (with more model cells) is used, we would
expect on the average the model resolution for a cell at the same location will be reduced since
the cell size is smaller. Theoretically the sum of the elements in a column of the model
resolution matrix (equation 4.8) is equals to 1.0. The average value of the array elements in the
column would then be equals to 1.0/m, where m is the number of model cells. A more useful
measure to judge whether a model cell is resolved is the ratio the cell resolution value (i.e. the
diagonal element of the resolution matrix) to this average value. If a value of about 10 is used
for this ratio (model resolution index), then the average maximum depth of survey is probably
about 12 to 14 meters. It is slight deeper between the 72 and 120 meters marks and slightly
shallower between the 24 and 72 meters marks (Figure 65b). This is similar to that shown by
the DOI index section using the 0.1 contour as the cut-off value (Figure 65c) The model
resolution sections shows a more gradual change with depth (and also laterally) in the
resolution values, compared to the DOI index section that has localized regions with high DOI
index values.
One advantage of the model resolution section is that it avoids the localized regions
with high DOI index values (that is sometimes caused by local high resistivity regions, or noise
in the data). It has a theoretical appeal in that it is less empirical in nature compared to the DOI
index method. One disadvantage of the model resolution method is that it requires an inversion
of a matrix with a computer time that is proportional to n3 (where n is the number of model
cells). This makes it impractical for models with more than about 100000 cells that could be a
significant limit for very large 3-D models with present PC technology. Research is being
carried out to extend the practical limit to beyond about 100000 cells. One possible solution is
to use GPU technology (Loke et al., 2010b).
Figure 65. Comparison between the (a) model resolution, (b) model resolution index and (c)
DOI index sections for the Landfill data set.
Figure 66. The model resolution plots for the three streamer configurations.
The resolution plot for the dipole-dipole array streamer (Figure 66a) shows the
maximum depth of investigation is about 4.5 meters. For long survey lines that consists of
hundreds of measurements the end effects are negligible so the resolution values with almost
flat contours near the center of the section gives an estimate of the maximum depth with
significant resolution values. The performance of gradient array streamer (Figure 66b) is much
poorer. The resolution plot shows the maximum depth of investigation is about 3 meters that is
significantly less than the dipole-dipole array streamer. The maximum depth of investigation
of the Wenner-Schlumberger array streamer (Figure 66c) is slightly less than the gradient array
streamer, between 2.5 and 3 meters, and much less than the dipole-dipole array streamer.
Overall, the dipole-dipole streamer has the best depth penetration but it has the lowest signal
strength that has to be compensated for in a practical design. The results from the model
resolution plots are significantly different from the cumulative sensitivity plots (Figure 34), so
some care must be exercised in using the results from the sensitivity plots.
programs, actually carries out an optimization (i.e. not a direct one-to-one inversion in the sense
it must have only one solution) in that it tries to reduce the difference between the calculated
and measured apparent resistivity values. If there is infinite data and a perfect fit between the
calculated and measured values, how the data is measured should not have an effect on the
results. However with real, noisy and limited data, how the data is measured does have an
effect. A model with 5% RMS error in the fit between the measured and calculated apparent
resistivity values with one data set might not give the same model as a 5% RMS error with
another data set although both might be from the same place. For this reason, the dipole-dipole
array gives an inversion model with much better resolution than the pole-pole array although
in theory the dipole-dipole values can be extracted from the pole-pole values.
Figure 67. An example of 3-D effects on a 2-D survey. (a) Apparent resistivity pseudosections
(Wenner array) along lines at different y-locations over (b) a 3-D structure shown in the form
of horizontal slices.
Figure 68. The 2-D sensitivity sections for the pole-dipole array with a dipole length of 1 meter
and with (a) n=6, (b) n=12 and (c) n=18. Note that as the ‘n’ factor increases, the zone of high
positive sensitive values becomes increasingly concentrated in a shallower zone below the P1-
P2 dipole.
(j). Increasing the electrode separation does not always increase the survey depth. It is
generally assumed that as the separation between the electrodes is increased, the region of the
subsurface that is ‘sensed’ by the array also increases. While this is true of most arrays, there
are certain important exceptions. In particular, this is not true of the pole-dipole array under
certain circumstances. In some surveys with the pole-dipole array, the separation between the
C1 current electrode and the P1-P2 dipole is increased in an effort to increase the depth of
survey by the array. However, if this is done with the P1-P2 dipole length (the ‘a’ factor in
Figure 4) kept at a constant spacing, certain interesting effects come into play. In section 2.5.8
this practice was strongly discouraged on the basis that the potential will decrease with the
square of the ‘n’ factor. This problem can be overcome by a combination of using higher
currents and more sensitive receivers. However, the problem caused by the change in the array
sensitivity pattern as the ‘n’ factor is monotonically increased is usually not taken into account.
The change in the sensitivity pattern when the ‘n’ factor changes from 1 to 6 was shown earlier
in Figure 30. Figure 68 shows what happens when the ‘n’ factor jumps from 6 to 12 to 18.
Here, the dipole length is kept constant at 1 meter. When ‘n’ is equal to 6 there are reasonably
high sensitivity values to a depth of about 3 to 4 meters between the C1 current and the P1
potential electrode. When ‘n’ is increased to 12, the zone of high sensitivity values becomes
increasingly more concentrated below the P1-P2 dipole in an even shallower region. This
means that the array with ‘n’ equals to 12 is in fact less sensitive to deeper structures than the
array with ‘n’ equals to 6. This effect is even more pronounced when ‘n’ is increased to 18.
Thus increasing the separation between the current electrode and the potential dipole,
while keeping the dipole length fixed, does not increase the survey depth of the array. It, in
fact, effectively decreases the depth of the region sensed by the array!
Figure 69 shows the apparent resistivity anomaly due to a small near-surface high
resistivity block for the pole-dipole array for ‘n’ values of up to 28. Note that the amplitude of
the high resistivity anomaly due to the near-surface block increases with the ‘n’ value, i.e. the
array becomes increasingly more sensitive to the near-surface block as the separation between
the electrodes increase. In field surveys with the pole-dipole array where the ‘n’ factor is
monotonically increased in the belief that this increases the survey depth, the pseudosection is
frequently dominated by a series of parallel slanting high-amplitude anomalies due to near-
surface inhomogeneities. The anomalies due to the near-surface structures frequently mask the
anomalies due to deeper structures that are of interest.
Figure 69. Example of apparent resistivity pseudosection with pole-dipole array with large ‘n’
values. Note that the anomaly due to a small near-surface high resistivity block becomes greater
as the ‘n’ factor increases. This means that the sensitivity of the array to the near-surface region
between the P1-P2 potential dipole becomes greater as the ‘n’ factor increases.
To increase the depth of penetration with the pole-dipole array, the ‘a’ dipole length
should be increased when the ‘n’ factor exceeds 6 to 8. This method was discussed in detail in
section 2.5.9. For example, instead of fixing the ‘a’ dipole length to 1 meter and increasing the
‘n’ factor to 28 in Figure 69, a more prudent approach is increase the ‘a’ spacing from 1 to 4
meters while ensuring the ‘n’ factor does not exceed 8.
Figure 70. Diagrammatic illustration of differences in objective function shapes for the pole-
pole array and dipole-dipole array data sets leading to different models obtained from
optimization routine.
5 I.P. inversion
5.1 Introduction
One of the more recent developments in the instrumentation for electrical imaging
surveys has been the addition of Induced Polarization (I.P.) capability in the multi-electrode
resistivity meter system. Many of the early 2-D surveys were resistivity and I.P. surveys carried
out using conventional 4-electrode systems in the 1950's onwards for mineral exploration,
particularly for conductive sulfide ore bodies. Quantitative interpretation of such historical data
was rather limited due to the limited computing facilities available at that time. Such historical
data provides an interesting source of data for testing modern 2-D and 3-D inversion software.
Re-interpretation of such old data to produce quantitative models sometimes has shed new light
on the geological structures.
One of the distinctive characteristics of the I.P. method has been the different
parameters in the time and frequency domains used to represent the I.P. effect. The following
section briefly discusses the I.P. phenomena and the different I.P. parameters. This is followed
by a few exercises in the inversion of IP data with the RES2DINV program. The data format
used by the RES2DINV and RES3DINV programs is described in their respective manuals. As
such, we will not cover it here.
Figure 71. The I.P. values for some rocks and minerals.
Figure 72. The Cole-Cole model. (a) Simplified electrical analogue circuit model (after Pelton
et al. 1978). = resistivity, m = chargeability, = time constant, c = relaxation constant. (b)
Amplitude and phase response to sine wave excitation (frequency domain). (c) Transient
response to square wave current pulse (time domain). Most I.P. receivers measure the integral
of the decay voltage signal over a fixed interval, mt, as a measure of the I.P. effect.
One mathematical model that attempts to explain variation of resistivity with frequency
observed in the IP method is the Cole-Cole mode (Pelton et al., 1978), which is defined by
1
( ) 0 1 m1
c (5.1)
1 i
where 0 is the DC resistivity, m is the chargeability, is the angular frequency (2f), is a
time constant and c is the exponent or relaxation constant. While the DC resistivity and
chargeability determine the behavior of the material at very low and very high frequencies, the
variation of the amplitude and phase curves at intermediate frequencies are also affected by the
time and relaxation constants.
The time constant factor has a large range, from 0.01 second to several thousand
seconds. The relaxation constant factor c is bounded by 0.0 to 1.0, with values frequently
between 0.2 and 0.7. Much of the earlier work was on the use of spectral I.P. (SIP)
measurements to differentiate between different types of conductive minerals for mining
purposes (Van Voorhis et al.,1973; Zonge and Wynn, 1975; Pelton et al., 1978; Vanhalla and
Peltoniemi, 1992). More recently, attempts have been made to use the SIP method for
environmental surveys, such as in the detection contaminants (Vanhala et al., 1992). Figure 72
shows a simplified electrical analogue circuit for the Cole-Cole model, together with typical
response curves in the frequency and time domain.
mt 1870
0.15
,
Vs dt
(5.2)
VDC
where the integration is carried out from 0.15 to 1.1 seconds after the current cut-off. The
chargeability value is given in milliseconds. The chargeability value obtained by this method
is calibrated (Summer 1976) so that the chargeability value in msec. has the same numerical
value as the chargeability given in mV/V. In theory, the chargeability in mV/V has a maximum
possible value of 1000.
I.P. surveys have traditionally been used in the mineral exploration industry,
particularly for metal sulfides, where heavy electrical generators producing high currents of the
order of 10 Amperes are used. The apparent I.P. values from such surveys are usually less than
100 msec. (or mV/V). One recent development is the addition of I.P. capability to battery based
systems used in engineering and environmental surveys where currents of 1 Ampere or less are
normally used. An accompanying phenomenon is the observation of I.P. values of over 1000
msec. (or less than -1000 msec.) in some data sets. Such values are almost certainly caused by
noise due to a very weak IP signal. To check whether such high I.P. values are real, first check
the apparent resistivity pseudosection. If it shows unusually high and low values that vary in
an erratic manner, the data is noisy. If the apparent resistivity values are noisy, then the apparent
I.P. values are almost certainly unreliable. Next check the apparent I.P. pseudosection. If the
apparent I.P. values show an erratic pattern (frequently with anomalous values lined up
diagonally with an apex at a doubtful electrode), then the I.P. values are too noisy to be
interpretable. There has been some recent work on improving the reliability of I.P.
measurements made with the multi-electrode type of systems (Dahlin and Leroux, 2012).
(b) Frequency domain I.P. measurements
In the frequency-domain methods, the apparent resistivity is measured at least two
frequencies, for example at 1 Hz. and 10 Hz. The higher frequency is usually set at 10 times
time the lower frequency. One commonly used frequency domain IP unit is the percent
frequency effect, PFE, which is defined by
L H
PFE .100 (5.3)
H
Another closely related unit that is also commonly used is the metal factor MF which is defined
by
L H
MF .2 .10 5 (5.4)
L . H
Another common frequency domain IP unit is the phase angle, . It is the phase shift between
the transmitter current and the measured voltage, and the unit commonly used is milliradians.
This has the advantage that only measurement at a single frequency is need, but the current
circuit must be coupled with the potential measuring circuit in some way so that the phase shift
between the measured potential and the input current can be determined. This is not a problem
in the multi-electrode type of system, such as the Pasi Polares system, where both circuits are
in the same unit.
(c) Relationship between the time and frequency domain IP units
From the measurement of the amplitude and phase spectrum of porphyry copper
mineralization, Van Voorhis et al. (1973) proposed the following equation to describe the
observed spectra.
K j b (5.5)
K is constant and b is a measure of the IP effect. It is a positive number of more than 0 and less
than 0.1. This is also known as the constant phase model (Weller et al., 1996). By using the
above model, the following relationship between the different IP units and the b parameter
were derived.
PFE 10b 1.100
= 1571 b (5.6)
mt 1320b
These relationships provide a numerical link between the different IP units (Van Voorhis et al.,
1973; Nelson and Van Voorhis, 1973). An alternative relationship is given in the paper by
Kemna et al. (1997).
There have been a large number of 2-D I.P. surveys published over the years, so there
is certainly no lack of data to test the IP inversion with the RES2DINV program. Most of the
newer multi-electrode systems now come with an IP option. However, the data from many of
these systems is rather noisy due to the limited current that can be injected into the ground. For
reasonable IP data quality, a current of at least 0.5 Amperes is probably necessary. For
environmental and engineering surveys, the most useful application of the IP data is probably
in differentiating between sand and clay sediments.
Figure 73 below shows the inversion results from both I.P. inversion methods for the
data from the Magusi River ore body (Edwards, 1977) survey. The models obtained by both
methods are generally similar. The ore body shows up as a low resistivity body with high IP
values near the middle of the survey line in the model sections. The blocky inversion method
was used to sharpen the boundary between the ore body and the surrounding country rocks in
the model. The massive sulfide ore body shows as a prominent low resistivity and high I.P.
structure. The complex resistivity method gives a model with slightly higher I.P. values.
Figure 73. Magusi River massive sulphide ore body inversion models. (a) Apparent resistivity
pseudosection. Resistivity inversion models obtained using the (b) perturbation and (c)
complex resistivity methods. (d) Apparent I.P. (metal factor) pseudosection. I.P. models
obtained using the (e) perturbation and (f) complex resistivity methods.
Figure 74. Sketch of separated cable spreads setup used (after Dahlin and Loke, 2015).
Figure 75. Resistivity and chargeability pseudosection from field demo at 3rd IP workshop at
Ile d’Oleron (after Dahlin and Loke, 2015).
6 Cross-borehole imaging
6.1 Introduction
One of the most severe limitations of 2-D imaging surveys carried out along the ground
surface is the reduction in the resolution with depth. This is a fundamental physical limitation
that no amount of reconfiguration of the surface arrays or computer modeling can overcome.
In theory, the only way to improve the resolution at depth is to place the sensors (i.e. the
electrodes) closer to the structures of interest. This is not always possible, but when such
boreholes are present, cross-borehole surveys can give more accurate results than is possible
with surface surveys alone.
That has been many fine publications on such surveys; such as by Zhao et al. (1986),
Daily and Owen (1991), Sasaki (1992), LaBrecque et al. (1996), Slater et al. (1997, 2000),
Zhou and Greenhalgh (1997, 2000) and Wilkinson et al. (2006a, 2008). In the following
section, I will attempt to summarize the main results with regards to the choice of array
configurations for cross-borehole surveys. However, please refer to the original papers for the
details.
electrodes are below the surface of a half-space, the geometric factor is given by
rr '
k 4
r r'
where r’ is the distance of the reflected image of the current (Figure 77) from the potential
electrode.
Figure 76. The possible arrangements of the electrodes for the pole-pole array in the cross-
borehole survey and the 2-D sensitivity sections. The locations of the two boreholes are shown
by the vertical black lines.
Figure 77. A schematic diagram of two electrodes below the surface. The potential measured
at P can be considered as the sum of the contribution from the current source C and its image
C’ above the ground surface.
positive values.
Overall, many authors have positive remarks about this array. It provides better
resolution and is less sensitive to telluric noise (since the two potential electrodes are kept
within the survey area) compared to the pole-pole array. While in theory the resolution of the
array is slightly poorer than the bipole-bipole array, the potential values measured are
significantly higher.
Figure 78. The 2-D sensitivity patterns for various arrangements with the pole-bipole array.
The arrangement with (a) C1 and P1 in first borehole and P2 in second borehole, (b) C1 in the
first borehole and both P1 and P2 in the second borehole, (c) all three electrodes in the first
borehole and (d) the current electrode on the ground surface.
Figure 79. The 2-D sensitivity patterns for various arrangements of the bipole-bipole array. (a)
C1 and P1 are in the first borehole, and C2 and P2 are in second borehole. (b) C1 and C2 are
in the first borehole, and P1 and P2 are in second borehole. In both cases, the distance between
the electrodes in the same borehole is equal to the separation between the boreholes. The
arrangements in (c) and (d) are similar to (a) and (b) except that the distance between the
electrodes in the same borehole is half the spacing between the boreholes.
BORELUND.DAT – A pole- (1). Read in the file, and try an inversion with the default
pole array field data set from inversion parameters.
Lund University, Sweden. The (2). Next you can try different settings, such as the robust
survey was conducted to map inversion option.
fractures in a limestone-marl (3). The program also has an option to reduce the size of
formation (Dahlin pers. the model cells by half, which you might like to try.
comm.).
BORELANC.DAT – A bipole- (1). Read in the file, and run the inversion. The path of the
bipole array field data set from saline tracer is represented by regions with low resistivity
Lancaster University, U.K values. Can you identify the tracer in the model sections?
(Slater et al. 1997). The survey
was conducted to map the flow
of a saline tracer from the
ground surface through the
unsaturated zone.
Figure 80. Possible measurement sequences using the bipole-bipole array. Other possible
measurements sequences are described in the paper by Zhou and Greenhalgh (1997).
Figure 81. Several possible bipole-bipole configurations with a single borehole. (a) The C1
and C2 electrodes at depths of 3 and 2 meters respectively below the 0 meter mark. (b) The C1
and P1 electrodes at depths of 3 and 2 meters respectively below the 0 meter mark. (c). The C1
electrode is at a depth of 3 meters below the 0 meter mark while the C2 electrode is on the
surface. (d). The C1 electrode is at a depth of 3 meters below the 0 meter mark while the P1
electrode is on the surface.
Figure 82. A pole-bipole survey with a single borehole. The C1 electrode is at a depth of 3
meters below the 0 meter mark.
a conventional vertical borehole. This system is useful when it is necessary to image a region
of a long horizontal extent over a limited depth range. In a typical survey, the length of the line
is about 5 times the depth of the subsurface electrodes.
Figure 83. Test of optimized cross-borehole arrays with a synthetic model. (a) Two-layer test
model with conductive and resistive anomalies. Inversion models for (b) optimized data set
with all arrays, (c) 'standard' data set and (d) the reduced optimized data set that excludes arrays
with both current (or potential) electrodes in the same borehole. All the data sets have 1875
data points. The outlines of the rectangular blocks showing their true positions are also shown.
A method to generate optimized arrays for this type of survey configuration is described in
Loke et al. (2015b). An example of the improvement that can be obtained compared to
‘standard’ arrays generated manually using heuristic rules (Harro and Kruse, 2013) is shown
in Figure 85. The survey was carried out at the Geopark research site on University of South
Florida campus in west-central Florida, United States. The site is characterized by karstified
limestone bedrock overlain by about 5 meters of overburden soils consisting of granular sands
over more cohesive sandy clay and clay soil with clay content generally increasing with depth
(Loke et al., 2015b). Depths to contacts were available for standard penetration tests (SPTs),
cone penetration tests (CPTs) and GPR data (Stewart and Parker, 1992). A deep array of 14
electrodes was implanted at 7.62 m below ground surface with an electrode spacing of 4 m
with a matching set of 14 electrodes on the surface directly above the implanted electrodes.
The models for the optimized data sets (Figure 85b and c) shows better agreement with the
known geology from the geotechnical and GPR data, in terms of the lower boundary of the top
sandy layer and the positions of the cavities, compared to the ‘standard’ arrays model (Figure
85a).
Figure 84. Schematic diagram of the MERIT method with the electrodes are planted along the
surface and directly below using the direct push technology (after Harro and Kruse, 2013).
Figure 85. Inversion models for the different data sets for the data collected with electrodes at
surface and 7.62 m depth with 4 m horizontal spacing. Models for the (a) standard arrays (405
data points), optimized arrays with (b) 403 and (c) 514 data points.
7.1 Introduction
Here we will look at a number of examples from various parts of the world to give
you an idea of the range of practical survey problems in which the electrical imaging method
has been successfully used.
Figure 86. Landslide field example, Malaysia. (a) The apparent resistivity pseudosection for
a survey across a landslide in Cangkat Jering and (b) the interpretation model for the
subsurface.
Due to the nature of such sites, the subsurface is often very complex and is a challenging
target for most geophysical methods. The survey for this example was carried out on a
derelict industrial site where leachate was known, from a small number of exploratory wells,
to be moving from a surface waste lagoon into the underlying sandstones (Barker, 1996).
Eventually the leachate was seen seeping into a nearby stream. However, the extent of the
subsurface contamination was not known.
An electrical imaging survey was carried out along an old railway bed between the
lagoon and the stream. The metal railway lines had been removed except for short lengths
embedded in asphalt below a large metal loading bay. In the apparent resistivity pseudosection
(Figure 87a), the area with contaminated ground water shows up as a low resistivity zone to
the right of the 140 meters mark. The metal loading bay causes a prominent inverted V shaped
low resistivity anomaly at about the 90 meters mark. In the inversion model (Figure 87b), the
computer program has managed to reconstruct the correct shape of the metal loading bay near
the ground surface. There is an area of low resistivity at the right half of the section that agrees
with what is known from wells about the occurrence of the contaminated ground water. The
plume is clearly defined with a sharp boundary at 140 meters along the profile. The
contaminated zone appears to extend to a depth of about 30 meters.
Figure 87. Industrial pollution example, U.K. (a) The apparent resistivity pseudosection from
a survey over a derelict industrial site, and the (b) computer model for the subsurface.
anomaly is detected below the 200 ft. mark, which is probably a hole in the lower clay layer
(Figure 88b). This feature falls in an area in the pseudosection where there is an apparent gap
in the data.
However, a plot of the sensitivity value of the cells used in the inversion model shows
that the model cells in the area of the high resistivity body have higher sensitivity values (i.e.
more reliable model resistivity values) than adjacent areas at the same depth with more data
points in the pseudosection plot (Figure 88c). This phenomenon is due to the shape of the
contours in the sensitivity function of the dipole-dipole array (Figure 24), where the areas with
the highest sensitivity values are beneath the C1C2 and P1P2 dipoles, and not at the plotting
point below the center of the array. This example illustrates the danger of only using the
distribution of the data points in the pseudosection to constrain the position of the model cells
(Barker, 1992; Loke and Barker, 1996a). If the model cells are placed only at the location of
the data points, the high resistivity body will be missing from the inversion model, and an
important subsurface feature would not be detected!
Figure 88. Mapping of holes in a clay layer, U.S.A. (a) Apparent resistivity pseudosection for
the survey to map holes in the lower clay layer. (b) Inversion model and (c) sensitivity values
of the model cells used by the inversion program.
Figure 89. Water infiltration mapping, U.K. (a) The apparent resistivity and (b) inversion
model sections from the survey conducted at the beginning of the Birmingham infiltration
study. This shows the results from the initial data set that forms the base model in the joint
inversion with the later time data sets. As a comparison, the model obtained from the inversion
of the data set collected after 10 hours of irrigation is shown in (c).
The water distribution is more easily determined by plotting the percentage change in
the subsurface resistivity of the inversion models for the data sets taken at different times
(Figure 90) when compared with the initial data set model. The data set collected at 5 hours
after the pumping began shows a reduction in the resistivity (of up to over 50 percent) near the
ground surface in the vicinity of the 24 meters mark. The near-surface low resistivity zone
reaches its maximum amplitude after about 10 hours when the pumping was stopped (Figure
90b). Twelve hours after the pumping was stopped, the low resistivity plume has spread
downwards and slightly outwards due to infiltration of the water through the unsaturated zone.
There is a decrease in the maximum percentage reduction in the resistivity values near the
surface due to migration of the water from the near surface zone. This effect of the spreading
of the plume becomes increasingly more pronounced after 24 hours (Figure 90d) and 36 hours
(Figure 90e) due to further migration of the water. Note that the bottom boundary of the zone
with approximately 20 percent reduction in the resistivity values tends to flatten out at a depth
of about 3 meters (Figure 90e) where the plume from the surface meets the water table.
Figure 90. Time-lapse sections from the infiltration study. The sections show the change in the
subsurface resistivity values with time obtained from the inversion of the data sets collected
during the irrigation and recovery phases of the study.
Figure 91. Hoveringham pumping test, U.K. (a) The apparent resistivity pseudosection at the
beginning of the test. The inversion model sections at the (b) beginning and (v) after 220
minutes of pumping.
Figure 92. Percentage relative change in the subsurface resistivity values for the Hoveringham
pumping test. To highlight the changes in the subsurface resistivity, the changes in the model
resistivity are shown. Note the increase in the model resistivity below the borehole with time.
Figure 92 clearly shows the increase in the zone with higher resistivity values with time due to
the extraction of the water. By using Archie’s Law, and assuming the water resistivity does not
change with time, we can estimate the change in the water saturation values. The decrease in
the water saturation level within the aquifer, or desaturation values, is shown in Figure 93. As
Archie’s Law assumes that the conduction is due to the water content alone, the desaturation
values are likely to be lower than the true values if there is significant clay content.
Figure 93. Use of Archie’s Law for the Hoveringham pumping test. Sections showing the
relative desaturation values obtained from the inversion models of the data sets collected during
the different stages of the Hoveringham pumping test. Archie’s Law probably gives a lower
limit for the actual change in the aquifer saturation.
To emphasize the boundary between the soil layer and the bedrock, the robust inversion
option was used (section 4.3). The inversion model is shown in Figure 94b. The thickness of
the lower resistivity weathered layer is generally about 10 to 20 meters. There is a narrow
vertical low resistivity zone with a width of less than 20 meters below the 190 meters mark that
is probably a fracture zone in the bedrock. A borehole well that was placed at the 175 meters
mark that lies just at the edge of the fracture zone. It had yields that were somewhat lower than
expected (Acworth, 1987). In such a situation, the 2D resistivity model would be useful to
pinpoint the exact location of the center of the fracture zone to improve the yield from the
borehole. The placement of the well was largely based on resistivity and EM profiling data,
and many years before 2D resistivity inversion software and fast microcomputers were widely
available.
Figure 94. Groundwater survey, Nigeria. (a). Apparent resistivity pseudosection. (b) The
inversion model with topography. Note the location of the borehole at the 175 meters mark.
As a final note, it is possible to invert data collected with the Wenner Alpha, Beta and
Gamma arrays along the same line simultaneously with the RES2DINV program as a single
data set. This can be done by using the "non-conventional array" option in the program where
the positions of all the four electrodes in an array are explicitly specified. This might be an
interesting method to combine the advantages of the different variations of the Wenner array.
the near surface lithology of the riverbed where there were plans to lay a cable. The data set
has a total of 7479 electrode positions and 6636 data points, whereas the inversion model used
has 19936 cells. On a 3.2 Ghz Pentium 4 computer, it took slightly less than 2 hours to process
this data set. On a newer PC, and using the appropriate software settings, it will probably takes
less than 15 minutes (see section 4.8)!
Figure 95. The inversion model after 4 iterations from an underwater riverbed survey by Sage
Engineering, Belgium.
In the inversion model (Figure 95b), most of the riverbed materials have a resistivity of
less than 120 m. There are several areas where the near-surface materials have significantly
higher resistivities of over 150 m. Unfortunately, geological information in this area is rather
limited. In the high resistivity areas, the divers faced problems in obtaining sediment samples.
The lower resistivity materials are possibly more coherent sediments (possibly sand with
silt/clay), whereas the higher resistivity areas might be coarser and less coherent materials.
Figure 96. Thames River (CT, USA) survey with floating electrodes. (a) The measured
apparent resistivity pseudosection. Inversion models obtained (b) without constraints on the
water layer, and (c) with a fixed water layer.
largest deposits are the Athabasca Oil Sands with over 200 Gm 3 of reserves. The resistivity log
across one of the oil sands bodies (Figure 97b) shows high resistivity values approaching 1000
ohm.m for the tar sands as compared to less than 30 ohm.m for the sediments, making the
resistivity method a potentially good exploration approach. An example of a resistivity section
obtained from a 2D survey (Figure 97c) shows prominent high resistivity zones associated with
the oil sands.
Figure 97. Survey to map oil sands, Alberta, Canada. (a) Location of major tar sands deposits
in Alberta, Canada. (b) Example of resistivity log and geologic column of Athabasca oil sands.
(c) 2-D resistivity model from imaging survey (Kellett and Bauman, 1999).
Besides these examples, 2-D imaging surveys have been carried for many other
purposes such as the detection of leakage of pollutants from landfill sites, areas with undulating
limestone bedrock, mapping of the overburden thickness over bedrock (Ritz et al., 1999),
leakage of water from dams, saline water intrusion in coastal aquifers, freshwater aquifers
(Dahlin and Owen, 1998), monitoring of groundwater tracers (Nyquist et al., 1999) and
mapping of unconsolidated sediments (Christensen and Sorensen, 1994). The resistivity
imaging method has also been used in underwater surveys in lakes and dams.
Figure 99. Two possible measurement sequences for a 3-D survey. The location of potential
electrodes corresponding to a single current electrode in the arrangement used by (a) a survey
to measure the complete data set and (b) a cross-diagonal survey.
In some cases, 3-D data sets are constructed from a number of parallel 2-D survey lines
(section 8.3). Ideally there should be a set of survey lines with measurements in the x-direction,
followed by another series of lines in the y-direction. The use of measurements in two
perpendicular directions helps to reduce any directional bias in the data.
However, in some cases, only the data from a series of survey lines in one direction is
available. This is particularly common if the surveys were originally conducted to provide 2-
D images. Sometimes the spacing between the “in-line” electrodes is significantly smaller than
the spacing between the lines. One important question is the maximum spacing between the
lines that can be used for the data to be still considered “3-D”. A useful guide is the 3-D
sensitivity plot. Figure 100 shows the sensitivity values on horizontal slices through the earth.
The electrodes are arranged along the 0 and 1 meter marks along the x-axis. Near the surface,
there is an approximately circular region with negative sensitivity values in the top two slices
at depths of 0.07 and 0.25 meter. The zone with the largest sensitivity (using the 4 units
sensitivity contour line as a guide) extends in the y-direction to slightly over half the electrode
spacing. This means to get a complete 3-D coverage, if the measurements are only made in the
x-direction, the spacing between the lines should not be much more than the smallest electrode
spacing used.
The papers by Bentley and Gharibi (2004) and Gharibi and Bentley (2005) give a fairly
comprehensive discussion on constructing 3-D data sets from 2-D survey lines.
Figure 100. 3-D sensitivity plots for the pole-pole array. The plots are in the form of horizontal
slices through the earth at different depths.
The pole-pole array has two main disadvantages. Firstly it has a much poorer resolution
compared to other arrays. Subsurface structures tend to be smeared out in the final inversion
model. The second disadvantage, particularly for large electrode spacings, is that the second
current electrode and potential electrode must be placed at a sufficiently large distance from
the survey grid. Both disadvantages have been discussed in detail in Section 2.5.7. Park and
Van (1991) who used this array for a field experiment found that about 15% of the
measurements did not satisfy reciprocity because the contributions from the remote electrodes
were significant. In general, this probably affects the readings with the larger spacings.
Figure 101. 3-D sensitivity plots for the pole-dipole array with n=1 in the form of horizontal
slices through the earth at different depths. The C1 electrode is the leftmost white cross.
Figure 102. 3-D sensitivity plots for the pole-dipole array with n=4 in the form of horizontal
slices through the earth at different depths.
Figure 103 and Figure 104 shows the sensitivity patterns for the dipole-dipole array
when the “n” factor is equal to 1 and 4 respectively. The sensitivity values have small but
negative values outside the immediate vicinity of the array. Another interesting feature is that
the sensitivity contours tend to be elongated in the y-direction, particularly for the larger “n”
value. The 4 unit sensitivity contour extends to about 0.6 times the array length in the y-
direction, or about 1.8 times the unit electrode spacing. This means that the array is more
sensitive to structures off the array axis compared to the pole-pole and pole-dipole arrays. This
feature is troublesome in 2-D surveys, but might be advantageous in 3-D surveys.
The off-axis elongation of the sensitivity contours agrees with the observation by
Dahlin and Loke (1997) that the dipole-dipole array is more sensitive to 3-D effects compared
to other common arrays. This factor is important when the dipole-dipole array is used in 2-D
imaging surveys where it is assumed that the subsurface geology is 2-D. In many cases, 3-D
data sets for the dipole-dipole arrays are constructed from a number of parallel 2-D survey
lines, particularly from previous surveys. Due to the elongated sensitivity pattern, the dipole-
dipole array can probably tolerated a larger spacing between the survey lines (to about three
times the inline unit electrode spacing) and still contain significant 3-D information.
In closing this section on four electrodes arrays, the sensitivity patterns for the Wenner
alpha (Figure 105), Schlumberger (Figure 106) and Wenner gamma arrays (Figure 107, the
Wenner beta is the dipole-dipole with a “n” value of 1 as shown in Figure 101) are shown. The
sensitivity contours for the Wenner alpha array, outside of the immediate vicinity of the
electrodes, are elongated in the direction of the line of electrodes. This means that the Wenner
alpha array is less sensitive to off-line structures than the dipole-dipole array, i.e. it is less
sensitive to 3-D. This agrees with empirical observations by Dahlin and Loke (1997). The
sensitivity pattern for the Wenner-Schlumberger array (Figure 106) is similar to that for the
Wenner alpha array except for a slight bulge near the center of the array. The sensitivity
patterns for the Wenner gamma array (Figure 107) show characteristic bulges near the C1 and
C2 electrodes that were observed earlier in the 2-D sensitivity sections (Figure 23). Thus it is
expected to be more sensitive to 3-D structures near the C1 and C2 electrodes.
Figure 103. 3-D sensitivity plots for the dipole-dipole array with n=1 in the form of
horizontal slices through the earth at different depths. The C2 electrode is the leftmost white
cross.
Figure 104. 3-D sensitivity plots for the dipole-dipole array with n=4 in the form of
horizontal slices through the earth at different depths. The C2 electrode is the leftmost white
cross.
Figure 105. The 3-D sensitivity plots for the Wenner alpha array at different depths.
Figure 106. The 3-D sensitivity plots for the Wenner-Schlumberger array with the n=4 at
different depths.
Figure 107. The 3-D sensitivity plots for the Wenner gamma array at different depths.
measurements in the x- and y- directions only, without the diagonal measurements. This is
particularly common if the survey is made with a system with a limited number of independent
electrodes, but a relatively large grid is needed. A roll-along procedure using only 3 parallel
cables is described in Dahlin and Bernstone (1997).
In some cases, measurements are made only in one direction. The 3-D data set consists
of a number of parallel 2-D lines. The data from each 2-D survey line is initially inverted
independently to give a series of 2-D cross-sections. The measured apparent resistivity values
from all the lines can also be combined into a 3-D data set and inverted with RES3DINV to
give a 3-D picture. While the quality of the 3-D model is expected to be poorer than that
produced with a complete 3-D survey, such a “poor man’s” 3-D data set could reveal major
resistivity variations across the survey lines (see also section 8.5). Until multi-channel
resistivity instruments are widely used, this might be the most cost-effective solution to extract
some 3-D information from 2-D surveys. This arrangement might be particularly useful for
surveys with the dipole-dipole array that can tolerate larger spacings between the survey lines.
Figure 108. Using the roll-along method to survey a 10 by 10 grid with a multi-electrode
system with 50 nodes. (a) Surveys using a 10 by 5 grid with the lines orientated in the x-
direction. (b) Surveys with the lines orientated in the y-direction.
BLOCK15.MOD – A 15 by 15 (1). Read in the file, and then run the “Calculate” step to
grid model with several calculate the potential values.
rectangular prisms. (2). Next choose the “Display apparent resistivity” option
under the “Edit/Display” menu. Take a look at the
pseudosections for a few arrays.
(3). Try using the “Edit model” option to change the model,
and then recalculate the potential values. Check out the
effect of your changes on the apparent resistivity
pseudosections.
Figure 109a show an example of a 3-D model with a 15 by 15 survey grid (i.e. 255
electrodes). The model, which consists of four rectangular prisms embedded in a medium with
a resistivity of 50 m, is shown in the form of horizontal slices through the earth. The apparent
resistivity values for the pole-pole array (with the electrodes aligned in the x-direction) are
shown in the form of horizontal pseudosections in Figure 109b. Note the low resistivity block
with a resistivity of 10 m near the centre of the grid that extends from a depth of 1.0 to 3.2
meters. For measurements with the shorter electrode spacings of less than 4 meters this block
causes a low resistivity anomaly. However, for electrode spacings of greater than 6 meters, this
low resistivity prism causes a high resistivity anomaly! This is an example of “anomaly
inversion” which is caused by the near-surface zone of negative sensitivity values between the
C1 and P1 electrodes (Figure 100).
Figure 109. A 3-D model with 4 rectangular prisms in a 15 by 15 survey grid. (a) The finite-
difference grid. (b) Horizontal apparent resistivity pseudosections for the pole-pole array with
the electrodes aligned in the x-direction.
ROOT7_REMOTE.DAT – The Try inverting this data file and see whether the remote
same data set with the location of electrodes have a significant effect of the results.
the C2 and P2 electrodes included.
BLOCKS11T.DAT – A synthetic (1). Read in the file, and then invert the data set.
data set with topography. It has a 11 (2). Invert again using the “Robust inversion” option.
by 11 grid.
Figure 111. Inversion models for the Vetlanda landfill survey data set. (a) Using standard
inversion settings. (b) With a higher damping factor for the topmost layer. (c) Using diagonal
roughness filter in the horizontal (x-y) directions. (d) Using diagonal roughness filters in the
vertical (x-z and y-z) directions as well. (e) Using the roughness filter in all directions.
Figure 112. Types of 3-D roughness filters. (a) With components in the x- and y- directions
only for the horizontal filter. (b) With components in the diagonal directions in the x-y plane
for the horizontal filter. (c) Applying the roughness filter with the corner model cells as well.
Only two (out of eight) corner cells are shown.
In the final model, the diagonal components for the roughness filter are also applied
between the model cell and the cells adjacent to the corners of the cell (Figure 112e). Except
for the cells at the surface, bottom and sides of the model grid, there are eight such corner
neighboring cells surrounding each cell. The roughness filter is now applied between each
interior model cell and all the 26 cells surrounding it. The resulting inversion model shown is
in Figure 111e. In theory, this should reduce any bias in shape of the structures from the
alignment of model cells grid (Farquharson, 2008; Loke and Dahlin, 2010).
The model in Figure 111 has model cells with widths of 1 meter in the x direction and
2 meters in the y direction, i.e. elongated in the y direction. This is because the model
discretization follows the spacing between the electrode positions. It is possible the banding in
the y direction is partly caused by the model cells being elongated in the y direction. A possible
alternative method to reduce the banding is to use model cells that are of the same widths in
both directions, i.e. a selective model refinement in the y direction. This doubles the number of
model cells as well as the number of nodes in the finite-difference grid. Figure 113 shows the
inversion models using cells of 1 meter length in both the x and y directions. The inversion
model with the standard inversion settings (Figure 113a) does not exhibit the prominent
banding in the y direction (Figure 111a). Instead there is a slight banding in the x direction in
the top two layers with clustering of very low or high resistivity values along the survey lines.
This feature is removed when a higher damping factor is used for the first layer (Figure 113b).
The boundary of the low resistivity zone (marked by the light blue color) between x distances
of 15 to 27 meters in the top three layers in this model shows a 'stepped' pattern aligned in the
x or y direction. This is probably due to the use of a roughness filter with components in these
directions (Figure 112b). This feature is removed when the roughness filter with diagonal
components is used (Figure 112d) that does not bias the structure in the x and y directions
(Figure 113c). The inversion model with the diagonal filter components in the z direction as
well (Figure 113d) is also free of the 'stepped' boundary pattern.
Figure 113. Inversion models for the Vetlanda landfill survey data set using model cells of
equal lengths in the x- and y- directions. (a) Using standard inversion settings. (b) With a higher
damping factor for the topmost layer. (c) Using diagonal roughness filter in the horizontal (x-
y) directions. (d) Using diagonal roughness filters in the vertical (x-z and y-z) directions as well.
The resistivity of the metal well casing is much lower than the surrounding medium, so the
potential on the surface of the casing is essentially constant. The use of such ‘long’ electrodes
has been the subject of several papers (Schenkel and Morrison, 1994; Singer and Strack, 1998;
Ramirez et al., 2003; Daily et al., 2004). Many of these studies involve deep surveys where the
measurements are made using the metal cased wells as electrodes. One situation of interest for
environmental and engineering surveys is the use of both point and long electrodes in a survey
where existing metal wells can be used as the long electrodes (Rucker at al., 2007, 2010). It is
well known that well-to-well measurements involving vertical ‘long’ electrodes only have poor
vertical resolution. However, measurements between a ‘point’ and a ‘long’ electrode have
better vertical resolution. It might even help when the long electrode extends below the target
of interest.
Figure 114. Synthetic model for long electrodes survey. 3-D model using cells of low
resistivity (0.01 .m) that are marked in red to simulate cased wells.
Figure 114 shows a model with a low resistivity block of 20 .m at a depth of 2.9 to
6.0 meters. The maximum length of the block is 3 meter that is less than its average depth, thus
it will be a difficult target for normal surveys with point electrodes on the surface. Figure 115a
shows the inversion model for a pole-pole survey data set using the cross-diagonal
measurement sequence (Figure 99) using only point electrodes on the surface. Another data set
was generated for the same model with three of the point electrodes converted into long
electrodes. This was done by assigning a very low resistivity value to the cells around the point
electrode, as marked by the red cells in Figure 114. Two of the long electrodes extend below
the low resistivity block. The inversion model for the data set using the combination of point
and long electrodes (Figure 115b) shows slightly better resolution for the low resistivity block.
The lowest resistivity achieved in model is about 77 .m (Figure 115b) which is lower than
the value of 83 .m achieved by the model for the data set using only surface point electrodes
(Figure 115a).
Figure 115. Comparison of inversion models using point electrodes with and without long
electrodes. (a) Inversion model for pole-pole data set using only surface point electrodes. (b)
Inversion model for pole-pole data set using 121 point and 3 long electrodes.
8.9.2 Model grid optimization for large surveys with arbitrary electrodes positions
The data format that separates the model grid from the survey electrodes positions
(Figure 116d) has made it possible to model data from very complex survey layouts. This
flexibility comes at the price of higher computer memory requirements. Thus this option is
only available in the 64-bit RES3DINVx64 program. In this section, we will take a closer look
at the finite-difference and finite-element mesh used and ways to optimize it so as to reduce
the memory required.
Firstly, the program subdivides each model cell by four mesh lines in both the x and y
directions (Figure 117). When the electrode falls on a node location (at the intersection of the
x and y mesh lines), it can be directly modeled by that node. There are two alternatives to model
the effect of an electrode at an arbitrary position with the finite-difference and finite-element
methods when it does not fall on a mesh node. The first is to calculate the potential at the
electrode by interpolating the potentials (Figure 117a) at the four nearest nodes in the mesh
(and similarly replace a single current electrode node by four equivalent current sources). The
second method moves the nearest node to the location of the electrode using a distorted finite-
element mesh (Figure 117b). For both methods, it is important that two electrodes used in the
same array are sufficiently far apart in the mesh setup used. Figure 118a shows a problem when
a current electrode is close to a potential electrode in the same array when the interpolation
method is used. Both the current and potential electrodes will share a common node if they are
too close together. To avoid this problem, there should be at least two mesh lines between the
two electrodes so that they do not share a common node (Figure 118b).
It is possible to avoid the situation in Figure 118a by using smaller model cells.
However, using small model cells will increase the number of nodes in the finite-difference or
finite-element mesh. It will also increase the total number of model cells in the inversion model.
Both of this will increase the computer memory and time required for the inversion of the data
set. It might be necessary to make some adjustments to the model grid specified to obtain the
optimum balance.
Figure 117. Methods to model the effect of an electrode using the finite-difference and finite-
element methods.
Figure 118. The use of an appropriate mesh spacing to obtain sufficient accuracy for electrodes
in the same array that close together.
To illustrate this, we used the data from a survey at the Hanford site in Washington state, USA
where the waste material was stored in trenches and concrete cribs (Rucker et al., 2007).
Different resistivity survey phases were carried out using 2-D lines (Figure 119). The
distribution of the electrodes does not fit into a simple rectangular grid. The first attempt to
model the data uses a model grid with cell widths of 5 meters. This produced a model grid with
125x98 lines in the x-y directions. The inversion model had 16 layers that gave a total of 192448
model cells. While the inversion of the data set could be carried out on a PC with 24 GB RAM,
there were a number of electrode that were too close in the mesh, as marked by the light blue
points in the grid display (Figure 120a). The next model grid has cell widths of 4 meters (Figure
8.23b). This eliminated electrodes that were less than 2 mesh lines apart, but this produced a
model grid with 156x122 lines in the x-y directions. The inversion model has 300080 cells, and
the finite-difference mesh had more than 1.5 million nodes and the PC with 24 GB RAM was
not able to run the inversion. The third model grid with cell widths of 4 meters in the x-direction
where the electrodes were close together (on the right side of the mesh) and 5 meters elsewhere.
This gave a model grid with 136x98 lines. The model has 209520 cells that could be handled
by the PC with 24 GB RAM. The final inversion model obtained using this model grid setup
is shown in Figure 121. The linear features in the second to fifth layers are due to the trenches
and concrete cribs. The leakage zones are marked by the prominent low resistivity zones in the
fifth and deeper layers. Note that the bottom left corner of all the layers show relatively less
resistivity variation compared to the other sections. This is basically due to the fact it has no
data coverage as the survey lines do not cross this section.
Figure 119. Map with survey lines and infrastructure at the Hanford site.
Figure 120. Types of model grids for the Hanford survey data set. (a) Using a 5 meters spacing
model for the entire area. (b) Using a 4 meters spacing model for the entire area. (c) Using a
mixture of 4 and 5 meters spacing model grid.
(Loke et al., 2003) in the inversion of this data. The L-curve method was used to estimate the
optimum inversion damping factor. Figure 123 shows the inversion models for the
conventional arrays (with 3152 data points) and the optimized data (3154 data points) sets. The
optimized data set has slightly more points to maintain symmetry in the array configurations
used. The top two blocks are poorly resolved by the conventional arrays data set (Figure 123a)
compared to the optimized data set (Figure 123b). The width of the topmost block anomaly is
twice the actual size while its maximum resistivity is less than 240 ohm.m. In comparison, the
optimized arrays inversion model has a maximum resistivity of about 640 ohm.m and the
correct width. The optimized arrays also achieve significantly higher resistivity values at the
locations of the second and the third blocks. Note that although the data misfit for the optimized
data set (0.66%) is slightly higher than that for the conventional data set (0.34%) due to the
higher average geometric factor of the optimized arrays, the model resolution is better.
Figure 122. Arrangement of survey lines using a 3 cable system with the Abem SAS
instrument.
Figure 123. Inversions model for (a) combined Wenner-Schlumberger and dipole-dipole data
set, (b) optimized data set. The actual positions of the blocks are marked by black rectangles.
An interesting special case of a 3-D survey is when the electrodes are confined to the
perimeter of the survey area (Tejero-Andrade et al., 2015). The arrays used for such surveys
are commonly based on heuristic rules, and very often designed for a perimeter with sharp
corners such as a rectangle (Baker et al., 2001), and might not be applicable for perimeters with
smooth shapes such as a circle. An algorithm to generate optimized arrays for perimeters of
any shape is described in Loke et al. (2015d). As an example, for a survey with 40 electrodes
arranged in a square perimeter, the comprehensive data set has nearly 180 thousand arrays.
Figure 124a shows model resolution plots for the first 6 layers for the comprehensive data set
that shows the maximum possible resolution for this survey setup. The highest resolution
values are in the first layer near the electrodes and decrease rapidly with depth and towards the
centre of the survey region. The resolution sections for the ‘conventional’ arrays consisting of
the ‘L and Corner’ arrays (Tejero-Andrade et al., 2015) with 946 data points have much lower
resolution values, particularly in the deeper layers (Figure 124b). The optimized arrays
generated by the ‘Compare R’ method with the same number of data points have significantly
higher resolution values (particularly in the 5 th layer), although the resolution values are well
below that of the comprehensive data set (Figure 124c). Figure 125 shows model resolution
sections in a vertical plane across the middle of the area. It more clearly shows that the
optimized arrays have higher resolution values than the conventional arrays particularly at the
middle and at depth. This is clearly shown by the relative model resolution sections calculated
using the ratio of the resolution of the data set with that of the comprehensive data set (lower
row of sections in Figure 125). The ‘conventional’ arrays clearly has lower resolution values
near the center of the survey region and at depth. The larger optimized data set with 2000 data
points has an average resolution of 0.091 compared to 0.111 for the comprehensive data set
with almost 100 times more arrays.
Figure 124. Horizontal sections showing the model resolution for the (a) comprehensive data
set, (b) standard arrays and (c) optimized arrays. The electrode positions are marked by small
green crosses in the top layer in (a).
Figure 125. Vertical cross-sections showing the model resolution for the comprehensive,
‘standard’, small and large optimized data sets.
Figure 126. The synthetic test model with two 1000 ohm.m rectangular blocks (marked by
black rectangles) embedded in a 100 ohm.m background medium.
To illustrate the performance of the different arrays, a test model consisting of two high
resistivity rectangular blocks of 1000 ohm.m embedded within a 100 ohm.m medium ( Figure 126)
was used. The first test data set consists of the conventional 'L and Corner' arrays (Tejero-
Andrade et al., 2015) with the maximum geometric factor set at 4147 m that gives 946 data
points. Gaussian random noise (Zhou and Dahlin, 2003) with a mean amplitude of 1 milliohm
was added to the calculated resistance values before they were converted to apparent resistivity.
The 'L and Corner' arrays model (Figure 4a) detects the blocks with maximum resistivity values
of about 110 ohm.m but the shape of the deeper block is not clearly defined. The smaller
optimized data set model has a higher maximum value of about 119 ohm.m for the deeper block
(Figure 4b) with the maximum anomaly value that is closer to the true position in the upper
part of the block. However the top edge of the block is not well resolved with a small artifact
that extends to the surface. The poor vertical resolution is expected from the PSF plot (Figure
3e). The model for the larger optimized set with 2000 data points (Figure 4c) shows higher
values of 115 and 124 ohm.m at the positions of the top and bottom high resistivity blocks.
Although the optimized data sets have higher models misfits (1.5%) compared to the 'L and
Corner' arrays (0.5%) due to higher average geometric factors of the arrays used, the deeper
block is better resolved in the inversion models. The array optimization program took 616
seconds to generate the data set with 2000 arrays.
Despite being able to detect the larger block near the center of the survey grid, the
optimized data sets cannot clearly resolve the top of this structure. It was also observed by Tejero-
Andrade et al. (2015) that it was not possible to resolve the top of the structures used in their synthetic
models tests. Figure 128 shows plots of the point-spread-function (PSF) values (Oldenborger
and Routh, 2009) for a model cell at different depths located near the center of the square for
the comprehensive data set. All the PSF plots have a predominantly vertical pattern. The
regions with significant PSF values are spread out about 2 m laterally but about 4 m vertically.
As the comprehensive data set contains all the possible arrays, this indicates any survey with
the electrodes confined to the perimeter will be able to resolve horizontal boundaries better
than vertical boundaries. When the depth of the model cell is less than 2 m (or less than half
the distance from the nearest electrode) the highest PSF values occur near the surface. This is
the reason the anomalies in the inversion models reach the surface although the actual top
boundary of the structure is deeper. The maximum value of the PSF plot for the cell at 3.08 m
depth (Figure 128e) is located below the surface closer to the actual depth of the cell. Thus for
depths of more than half the distance from the nearest electrode the maximum model anomaly
value is expected to be below the surface although detecting the structure will more difficult
due to the loss of resolution with depth.
An example with a survey on a circular perimeter is next used to illustrate the general
nature of the array optimization algorithm. The array optimization algorithm has the advantage
over heuristic rules (based on a particular perimeter shape) as the same technique can be used
regardless of the shape of the perimeter. Figure 129 shows the model resolution sections for a
survey perimeter with 40 electrodes in a circle. As expected, the resolution decreases towards
the center of the survey area. The test model has two rectangular blocks of 10 ohm.m in a 100
ohm.m medium (Figure 130). The two blocks are detected in the inversion model for the
smaller optimized data set with 946 data points (Figure 130a). However, the top and bottom
boundaries of the deeper block are now well resolved. There is a slight improvement with the
larger optimized data set (Figure 130b). The left boundary of the top block and the bottom
boundary of the lower block are better resolved.
Figure 127. The inversion models for the synthetic model with two 1000 ohm.m rectangular
blocks (marked by black rectangles) embedded in a 100 ohm.m background medium.
Figure 128. Comprehensive data set point-spread-function plots on a vertical x-z plane located
at y=4.5 m for a model cell at different depths with centre at x=4.5m and y=4.5m.
Figure 129. Model resolution sections with circular perimeter for (a) comprehensive data set
with 180300 arrays, optimized data sets with (b) 946 and (c) 2000 arrays.
Figure 130. Inversion results for survey with a circular perimeter using optimized arrays with
(a) 946 and (b) 2000 data points.
In terms of use in actual field data inversion, the user defined model grid option
(probably with the arbitrary electrodes format) is likely to be used with the Res3dinvx64
program. This has the flexibility to adjust the x and y model grid spacings to fit the unit
electrode spacing used in the field survey. The model grid spacing should not be larger than
the unit electrode spacing used in the survey, and 4 nodes between adjacent model grid lines
should be used for the maximum accuracy. To effectively use more than 4 nodes between
adjacent electrode positions in the mesh, the user defined grid option allows using a grid
spacing that is smaller the unit electrode spacing. As an example, if a 50 meters electrode
spacing was used in the survey, and a model grid with 25 meter spacing is used, the finite-
element mesh will effectively have 8 nodes between adjacent electrode positions.
In designing the array configurations for a field survey, the geometric factor as well as
geometric factor relative error can be used to avoid arrays that are likely to be unstable. The
potential values measured by such unstable arrays are likely to be very small (and thus
overwhelmed by telluric noise), or very sensitive to errors in the electrode positions (this could
be a significant issue in field surveys in rugged and forested terrains), such that the data
measured is not meaningful and does not contribute significantly to resolving the subsurface
geology. This should be used to curb extreme field survey designs in an attempt to maximize
the data collection speed or coverage. In one data set I have come across, the arrays used have
geometric factors ranging from about 800 to 300 million m. Obviously measurements using
arrays with very large geometric factors are likely to be too noisy to be useful.
In selecting arrays to be used for field surveys, the following filtering steps are
suggested to eliminate arrays that are likely to be unstable.
1). Set an upper geometric factor limit, and remove arrays that exceed this limit. Arrays with
very large geometric factors will have very low potentials, and thus data from such arrays might
be too noisy to be useful. The upper geometric factor will depend on the characteristics of
equipment used (maximum current, signal processing features to improve signal-to-noise ratio
etc.) and the survey environment. A variation of this is to set a limit of the ratio of the maximum
geometric factor to the minimum geometric factor of the arrays in the data set to reflect the
dynamic range of signals that can be detected by the instrument.
2). Set an upper limit for the geometric factor relative error. This is to eliminate the offset type
of arrays that might have a reasonable geometric factors but likely to be unstable. A value of 4
to 10 can be used for most data sets.
In the Res3dinvx64 program, the above filtering methods are set using the ‘File-Cutoff
factor to filter data’ option. The following dialog box is shown on selecting this option.
Setting a value of 0.0001 for the geometric factor cutoff value allows the maximum geometric
factor to be 10,000 times the minimum geometric factor in the data points selected. This is
probably higher than the dynamic range of most resistivity-meter systems. For the geometric
factor relative error, a cutoff value of about 5 seems to work for most data sets.
3). Conduct a homogeneous earth model test using a finite-difference modeling program, such
as RES3DMODx64. In this case, calculate the apparent resistivity values for the proposed array
configurations for a homogeneous numerical model, such as with a resistivity of 100 ohm.m.
In theory, the calculated apparent resistivity values should be the same as the model resistivity
value. If the calculated apparent resistivity value differs by more than 20% from the model
value, the array is likely to be unstable. In any case, the forward modeling routine used in the
inversion cannot provide reasonably accurate apparent resistivity values for such arrays, so a
numerical interpretation is not viable using such data points.
In relation to this, some of the offset type of arrays used in 3-D surveys also use very
long distances between the current electrodes and potential dipole. This results in arrays with
very low potentials which are particularly sensitive to the boundary conditions used in the
finite-difference or finite-element forward modeling routine. In some cases, the small
inaccuracies in the boundary conditions can cause a net negative potential to be calculated
which in turn results in negative apparent resistivity values. One method to reduce the small
errors in the boundary conditions is to use more buffer nodes near the edges of the mesh used.
In the Res3dinvx64 program, this is selected using the ‘Change Settings – Forward Modeling
Settings – Change number of nodes’ option that brings up the dialog box shown below.
Using an extended boundary increases the number of nodes in the mesh, and thus increases the
computer time and memory required for the inversion of the data set. To satisfy these
requirements, a multi-core CPU and 64-bit system is required.
8.12 Not on firm foundations – inversion with shifting electrodes in 2-D and 3-D
It is implicitly assumed that the positions of the electrodes used in 2-D or 3-D surveys
are accurately known when interpreting the data set. In recent years, there has been significant
developments in geoelectrical monitoring surveys to detect temporal changes in the subsurface
(Chambers et al., 2014; Perrone et al., 2014), such as mapping of gas flows from landfills (Loke
et al., 2014a) and the monitoring of unstable slopes (Supper et al., 2014). The measurements
are repeated a number of times using the same set of electrodes. The electrode positions are
measured at the start of the campaign and possibly also at regular intervals. However, ground
movements sometimes occur between the times of electrode positions measurements. Thus the
electrode positions might not be accurately known for some data sets, making it necessary to
determine the shifts in the electrode positions from the resistivity data itself (Wilkinson et al.,
2015). A method was presented by Kim (2014) where both the subsurface resistivity and
electrodes positions are unknown variables to be determined by the least-squares optimization
method. The least-squares optimization equation (for a 2-D problem) is modified to the
following form.
G T
i
R d G i λ i V T R m V Δqi G Ti R d g i λ i V T R m V q i 1 (8.4)
r J W
where q x , G X , V αWx .
z Z βW
z
The combined Jacobian matrix G consists of J (the model resistivity sensitivity values), X and
Z (the sensitivity values for the changes in the x and z electrode positions) matrices. is a
damping factor vector and g is the data misfit vector. q is the model parameters vector
consisting of r (model resistivity) and x and z (electrode position) vectors. qi is the change
in the model parameters. W is the resistivity spatial roughness filter term. Wx and Wy are the
roughness filters for the electrode position vectors. Rd and Rm are weighting matrices used by
the L1-norm inversion method (Loke et al., 2003). and are the relative damping factor
weights for the shifts in the electrode positions.
Kim (2014) used the perturbation method (McGillivray and Oldenburg, 1990) to
calculate the spatial Jacobian matrices X and Z. For a survey line with e electrodes it will be
necessary to recalculate the potentials by resolving the finite-element capacitance matrix
equation 2e times. While this is possible for 2-D problems, it becomes impractical for 3-D
problems where the forward modeling routine is about two to three orders of magnitude slower.
Furthermore, it is necessary to calculate the partial derivatives in 3 directions (x, y, z). A much
faster method using a modification of the adjoint-equation approach (McGillivray and
Oldenburg, 1990) to calculate the X and Z Jacobian matrices is described in Loke et al. (2015c).
For 3-D problems, the adjoint-equation method is one to two orders of magnitude faster than
the perturbation method. For 3-D problems, it is probably three orders of magnitude faster. A
very large value of the spatial damping factors ( and ) of about 100 was also used by Kim
(2014) to reduce distortions in the electrode positions caused by near-surface resistivity
variations. However, it was found that using such large values for the and damping factors
tend to reduce the accuracy of the recovered electrode positions as well (Loke et al., 2015c).
A homogeneous half space is commonly used as the starting model for the optimization
algorithm. In a time-lapse survey the initial positions of the electrodes are usually accurately
measured, and can be treated as fixed parameters in the inversion of the initial data set. The
temporal changes in the subsurface resistivity are usually much less than the spatial variations
(Loke et al, 2014a). Thus the inversion model of the initial data set provides a good starting
model for the later time data sets. Loke et al. (2015c) modified the inversion algorithm to use
the model from a previous survey as the starting model which essentially eliminated distortions
in recovered positions of the x and z electrodes caused by near-surface resistivity variations.
Figure 131 shows a synthetic 2-D model used to illustrate the effect of shifts in the electrode
positions without the need to use large values for the spatial damping factors. The initial model
has 31 electrode positions with a uniform 1 m spacing on a flat surface (Figure 131a). Figure
131b shows the perturbed model with four changes. The electrode at the 5.0 m mark is shifted
0.3 m. to the right. The electrode at the 17.0 m mark is shifted 0.4 m. upwards. A 70 ohm.m
prism (depth of 1.0 m to the top) is added between the two existing prisms, and the 20 ohm.m
low resistivity prism is extended downwards by 0.73 m. The measured data consists of dipole-
dipole arrays with the 'a' dipole lengths ranging from 1 to 4 m., and the 'n' factor ranging from
1 to 6. This gives a total of 415 data points. Gaussian random noise (Zhou and Dahlin, 2003)
with a mean amplitude of 2.5 milliohm was added to the data before they were converted to
apparent resistivity values. The average noise level for the data set is about 1.0%. Figure 131b
shows the inversion models obtained assuming the electrodes are equally spaced on a flat
surface. The model for the perturbed data set has significant distortions (Figure 131b) below
the electrode at the 5m mark which has shifted horizontally, as well as below the 17 m mark
that has shifted vertically. Figure 131c shows the inversion model obtained when the electrodes
are allowed to shift. The model obtained from the inversion of the initial data set shown in
Figure 131a was used as the starting model. The distortions in the model have been removed
and the data misfit of 1.0% is close to the noise added. Note the upwards shift in the electrode
at the 17 m mark that matches shift in the synthetic model. In this inversion the spatial relative
damping factors and were set at 1.0.
Figure 132 shows a 3-D example with shifted electrodes. The initial model has the
electrodes arranged in a 29 by 13 rectangular grid with electrodes 1 meter apart on a flat surface.
The background medium has a resistivity of 100 ohm.m. There are two low (30 ohm.m) and
high (300 ohm.m) resistivity bands to demonstrate the effect of large near surface anomalies
on the inversion. There is a small near surface high resistivity rectangular prism (400 ohm.m)
and a deeper low resistivity prism (20 ohm.m). In the perturbed model, 3 of the electrodes are
shifted horizontally and 1 electrode vertically (Figure 132b). The resistivities of the small
rectangular prims are also changed to 25 and 350 ohm.m to simulate changes in the subsurface
resistivity. Voltage dependent random noise was added to the data which gave an average noise
level of about 1.1%.
Figure 131. The (a) initial and (b) perturbed synthetic test models with apparent resistivity
pseudosections and inversion models assuming a constant electrode spacing and flat surface.
(c) The inversion model obtained with the algorithm that allows the electrodes to shift.
Figure 133 shows the inversion results with fixed electrodes. The model for the initial
data set closely matches the true structure (Figure 133a) with a data misfit of 1.2% that is
slightly above the noise level added. The model for the perturbed data set show significant
distortions at the positions of the 4 shifted electrodes (Figure 133b) and a higher data misfit of
1.8%.
Figure 132. (a) 3-D synthetic test model with a rectangular survey grid. In the perturbed model,
the resistivities of the smaller blocks were changed from 400 and 20 ohm.m to 350 and 25
ohm.m. (b) The survey grid for the perturbed model. The four electrodes shifted are marked by
red circles. Electrode 1 was shifted 0.3 m in the x-direction, electrode 2 moved 0.3 m in the y-
direction, electrode 3 moved -0.2 m in both x and y-directions while electrode 4 was shifted
vertically upwards by 0.4 m.
Figure 133. Inversion models for the (a) initial and (b) perturbed data sets with fixed electrodes
in a rectangular grid. The true positions of the bands and prisms are marked by black lines. The
positions of the shifted electrodes are marked by small crosses in the top layer in (b).
Figure 134 shows the inversion models obtained when the electrodes positions are
allowed to change during the inversion using different values for electrodes relative damping
factor. The data misfit is similar to the added noise level with damping factors of 1.0 and 5.0
and the distortions in the resistivity structure seen in the model with fixed electrodes (Figure
134b) are removed. Using a high damping factor of 50 essentially fixes the electrode positions
(Figure 134c) giving a model similar to that with fixed electrodes (Figure 133b). However, the
recovered electrodes grid shows significant distortions with a damping factor of 1.0. It is
expanded outwards along the y-direction over the low resistivity band and compressed inwards
over the high resistivity band. This is probably because an increase in the spacing between the
electrodes reduces the measured resistance value in a similar way as a decrease in the
resistivity. The distortions are greatly reduced with a damping factor of 5.0. Figure 135a shows
the x-z surface profile along the y=10 m line that crosses electrode 4 (Figure 132b) that was
shifted upwards. The profile with a damping factor of 1.0 shows significant distortions that is
greatly reduce when it is increased to 2.5, and almost completely eliminated with a value of
5.0. Figure 136 shows the results obtained when the model obtained for the initial data set
(Figure 136a) is used as the starting model for the inversion of the perturbed data set. The
distortions with the lowest damping factor of 1.0 is greatly reduced (Figure 135b and Figure
136a) and almost eliminated with a damping value of 2.5 (Figure 135b). Using the initial data
set model essentially carries out the inversion using the change in the apparent resistivity values
which removes the effect of the common background resistivity structures.
Figure 134. Inversion models for the perturbed data set using different relative damping factors
for the electrodes positions vector with a homogenous half-space starting model.
Figure 135. Surface x-z profiles along the line y=10 m for inversions using a (a) homogenous
half-space and (b) initial data set starting models.
Figure 136. Inversion models for the perturbed data set using different relative damping factors
for the electrodes positions vector with the inversion model from the initial data set as the
starting model.
Figure 137. (a) Arrangement of electrodes in the Birmingham 3-D field survey. (b) Horizontal
and (c) vertical cross-sections of the model obtained from the inversion of the Birmingham
field survey data set. The locations of observed tree roots on the ground surface are also shown.
Figure 138 shows the use of the L-curve method (Farquharson and Oldenburg, 2004)
to estimate the optimum damping factor (section 1.4) in the inversion of this data set. The L-
curve method estimates the data misfit and model roughness for a large range of damping factor
values such as starting from one-hundredth to 100 times a trial value. The shape of the curve
obtained by plotting the model roughness versus the data misfit usually has an ‘L’ shape with
a distinct corner (Figure 138a). The damping factor value that corresponds to the maximum
curvature point of this curve is then chosen as the optimum damping factor value (Figure 138b).
Figure 138. Example L-curve plots. (a) A plot of the model roughness versus the data misfit
for the Birmingham survey data set for a number of damping factor values. (b) A plot of the
curvature of the curve in (a) for the different damping factor values. The ‘optimum’ damping
factor is at the maximum curvature point.
in Chambers et al (2006).
Figure 139. The 3-D model obtained from the inversion of the Lernacken Sludge deposit
survey data set displayed as horizontal slices through the earth.
Figure 140. A 3-D view of the model obtained from the inversion of the Lernacken Sludge
deposit survey data set displayed with the Slicer/Dicer program. A vertical exaggeration factor
of 2 is used in the display to highlight the sludge ponds. Note that the color contour intervals
are arranged in a logarithmic manner with respect to the resistivity.
Figure 141. Map of the Panama Canal region with the survey area marked (Noonan and
Rucker, 2011).
Figure 142. The model grid used for the Panama Canal floating electrodes survey with 20 by
20 meters cells. The location of the electrodes along the surveys lines are shown as colored
points.
Figure 143. The inversion model for the Panama Canal floating electrodes survey. The first 4
layers correspond to the water column.
Figure 144. Geological map of the Copper Hill area (White et al. 2001).
The 3-D I.P. model in Figure 146 shows two en-echelon north-south trends and two
approximately east-west trends forming an annular zone of high chargeability. The results from
existing drill-holes that had targeted the shallower part of the western zone agree well with the
resistivity and IP model. A drill-hole, CHRC58, intersected a 217m zone with 1.7 g/t gold and
0.72% copper coincided well an IP zone of greater than 35mV/V. The lower boundary of the
western zone with high chargeability coincides well with low assay results from existing drill-
holes. The eastern zone with high chargeability and resistivity values do not outcrop on the
surface and very little drilling has penetrated it. Further surveys, including drilling, is presently
being carried out.
Figure 145. Electrodes layout used for the 3-D survey of the Copper Hill area.
Figure 146. The I.P. model obtained from the inversion of the Copper Hill survey data set.
Yellow areas have chargeability values of greater than 35 mV/V, while red areas have
chargeability values of greater than 45 mV/V (White et al., 2001).
to other arrays (assuming these arrays used a sufficient long survey line to get the same
penetration depth). The pole-pole array offers field advantages of much shorter measurement
arrays and the ability to penetrate deeper than any other array for the same survey line length,
although there is the extra effort of maintaining two infinites (David Bingham, pers. comm.,
2009). The pole-dipole array has also been used for uranium exploration in the Athabasca basin,
and the 2-D and 3-D resistivity surveys played a significant role in the discovery of a very large
deposit (Donald Carriere, pers. comm., 2009) at Wheeler River near the existing McArthur
River deposit.
Figure 147. Location of uranium mines in the Athabasca basin, Saskatchewan (Bingham et al.,
2006).
Figure 149. Example of inversion model from the Midwest deposit area showing the resistivity
at a depth of about 200 meters (Bingham et al., 2006).
Figure 150. A complex time-lapse field survey example. (a) Map of Cripple Creek survey site.
(b) Overhead view of the inversion model grid with electrodes layout. (c) Iso-surface contours
for the -4% resistivity change at different times after the injection of the sodium cyanide
solution (that started at 2.8 hours from the first data set in snapshots used). t1= 1.1 hours, t2=
2.4 hours, t3= 3.7 hours, t4= 4.9 hours. (d) Overhead view of iso-surfaces.
Figure 151. Geological map of (a) south-west South Australia, (b) the Burra area, and (c) a
plot of survey electrodes and model cells layout.
The resistivity model (Figure 152a) shows a prominent north-south low resistivity
linear feature near the 1.8 km mark (x-axis) that corresponds to the Kingston Fault. The I.P.
anomaly (Figure 152b) in the northern part of the fault zone at depths of 100 to 200 meters
corresponds to the Eagle deposit prospect. The nature of the I.P. anomaly towards the bottom-
left edge of the deeper layers is uncertain as there is not much data coverage there. However it
lies in the Kingston Fault zone (Figure 151b) with reports of pyrites in a nearby bore.
Figure 153 show plots of the model resolution calculated using the resistivity and I.P.
Jacobian matrices, and the VOI using the model resistivity values. If a cutoff value of about 50
is used for the resistivity resolution index (Figure 153a), the maximum depth of investigation
is about 200 m. Not surprisingly, the highest resolution values are concentrated near the survey
lines, particularly around the group of shorter spacing lines in the northern third of the survey
area. This pattern is more pronounced in the I.P. resolution plots (Figure 153b). The I.P.
resolution sections have a shallower maximum depth of investigation than the resistivity
sections. This was confirmed by similar calculations for synthetic models with more uniform
data coverage. The VOI sections give a maximum depth of investigation of about 200 m. in the
southern half of the area below the longer survey lines (Figure 153c). The VOI plot show a
more complex pattern with local artefacts at several places and does not always increase
monotonically with depth (Figure 154c). The I.P. anomaly in the northern half of the survey
area lie in a region with higher resolution (and generally low VOI) values, so it is likely to be
real. The southern I.P. anomaly lie in a region with low resolution (and high VOI) values, so
its nature from the data alone is uncertain without independent confirmation.
Figure 152. Burra survey (a) resistivity and (b) I.P. inversion model layers.
Figure 153. The model (a) resistivity and (b) I.P. resolution index, and (c) VOI values. The red
arrows at the left side of (c) shows the position of the vertical slice shown in Figure 154.
Figure 154. Vertical cross-sections of the (a) resistivity model resolution index, (b) I.P
resolution index and (c) VOI in the X-Z plane 0.8 km north of the origin.
Some remarks about model reliability estimation - Some method of assessing the reliability
of the inversion results should be used. The sensitivity method (section 4.4) is the simplest but
it is difficult to pick a consistent cutoff value. The DOI/VOI method can be used for any data
set where it is possible to carry out an inversion. This plot normally has sharp boundaries that
make it easier to pick out the region with ‘reliable’ model values. A cut-off value of about 0.1
is usually used. However it is susceptible to local artefacts possibly caused by the local
optimization inversion method used. The model resolution method is more robust and less
affected by the inversion settings used. However it is computationally more demanding which
limits the size of the problem it can handle. Unlike the DOI/VOI plot, it shows a gradational
change in the resolution values that makes it more difficult to select a cut-off value.
Acknowledgments
Dr. Torleif Dahlin of Lund University in Sweden kindly provided the Odarslov Dyke,
Lernacken and Vetlanda data sets. The Grundfor data set was provided by Dr. Niels B.
Christensen of the University of Aarhus in Denmark and Dr. Torleif Dahlin. The Rathcroghan
data set was kindly provided by Dr. Kevin Barton and Dr. Colin Brown from data collected by
the Applied Geophysics Unit of University College Galway, Ireland. Many thanks to Richard
Cromwell and Rory Retzlaff of Golder Assoc. (Seattle) for the survey example to map holes in
a clay layer. Dr. Andrew Binley of Lancaster University kindly provided the interesting cross-
borehole data set. The Bauchi data set was provided by Dr. Ian Acworth of the University of
New South Wales, Australia. The Redas river survey data set was kindly provided by Jef
Bucknix of Sage Engineering, Belgium. Mr. Julian Scott and Dr. R.D. Barker of the School of
Earth Sciences, University of Birmingham provided the Clifton, the Tar Works and the 3-D
Birmingham survey data sets. Slicer/Dicer is a registered trademark of Visualogic Inc. I would
like to thank Paul Bauman and WorleyParsons Canada for permission to use the figures for the
Athabasca oil sands case history. I would like to thank David Bingham and Grant Nimeck of
Areva Resources Canada Inc. for permission to use the figures for the Athabasca uranium
deposit case history. The challenging Hanford, Panama Canal and Cripple Creek surveys data
set were provided by Dale Rucker of HydroGeophysics Inc. (USA). I would like to thank Kim
Frankcombe and Phoenix Copper for the Burra data set. Last, but not least, permission from
Arctan Services Pty. and Gold Cross Resources to use the Copper Hill 3-D example is
gratefully acknowledged.
References
Acworth, R.I., 1981. The evaluation of groundwater resouces in the crystalline basement of
Northen Nigeria. Ph.D. thesis, Univ. of Birmingham.
Acworth, R. I., 1987, The development of crystalline basement aquifers in a tropical
environment: Quarterly Journal Engineering Geology, 20, 265-272.
Alumbaugh, D.L. and Newman, G.A., 2000. Image appraisal for 2-D and 3-D electromagnetic
inversion. Geophysics, 65, 1455-1467.
Tejero-Andrade, A., Cifuentes-Nava, G., Chávez, R., Lopez-Gonzalez, A. and C. Delgado-Solorzano ,
2015. L- and CORNER-arrays for 3D electric resistivity tomography : an alternative for
geophysical surveys in urban zones. Near Surface Geophysics, 13, 335-367.
Aristodemou, E and Thomas-Betts, A, 2000. DC resistivity and induced polarisation
investigations at a waste disposal site and its environments. Journal of Applied
Geophysics, 44 (2-3), 275-302
Auken, E., and Christiansen, AV., 2004. Layered and laterally constrained 2D inversion of
resistivity data. Geophysics, 69, 752-761.
Baker, H.A., Djeddi, M., Boudjadja, A.G. and K. Benhamam, 2001. A different approach in
delineating near surface buried structures. EAGE 63rd Conference & Technical
Exhibition.
Barker, R.D., 1978. The offset system of electrical resistivity sounding and its use with a
multicore cable. Geophysical Prospecting, 29, 128-143.
Barker, R.D., 1989. Depth of investigation of collinear symmetrical four-electrode arrays.
Geophysics, 54, 1031-1037.
Barker R.D., 1992. A simple algorithm for electrical imaging of the subsurface. First Break 10,
53-62.
Barker R.D., 1996. The application of electrical tomography in groundwater contamination
studies. EAGE 58th Conference and Technical Exhibition Extended Abstracts, P082.
Barker, R. and Moore, J., 1998. The application of time-lapse electrical tomography in
groundwater studies. The Leading Edge, 17, 1454-1458.
Bauman, P., 2005. 2-D resistivity surveying for hydrocarbons – A primer. CSEG Recorder,
April, 25-33.
Bentley, L.R. and Gharibi, M., 2004. Two- and three-dimensional electrical resistivity
imaging at a heterogeneous remediation site. Geophysics, 69, 674-680.
Bernstone, C. and Dahlin, T., 1999. Assessment of two automated electrical resistivity data
acquisition systems for landfill location surveys : Two case histories. Journal of
Environmental and Engineering Geophysics, 4, 113-122.
Bingham, D., Nimeck, G., Wood, G. and Mathieson, T., 2006. 3D resistivity in the Athabasca
basin with the pole-pole array. 1 day workshop - Geophysical methods and techniques
applied to uranium exploration. SEG Annual General Meeting 2006, New Orleans.
Blome, M., Maurer, H. and Greenhalgh, S., 2011. Geoelectric experimental design — Efficient
acquisition and exploitation of complete pole-bipole data sets. Geophysics, 76, F15-F26.
Carpenter, E.W. and Habberjam, G.M., 1956. A tri-potential method of resistivity prospecting.
Geophysical Prospecting, 29, 128-143.
Chambers , J.E., Kuras, O., Meldrum, P.I., Ogilvy, R.O. and Hollands, J., 2006. Electrical
resistivity tomography applied to geologic, hydrogeologic, and engineering
investigations at a former waste-disposal site. Geophysics, 71, B231-B239.
Chambers J E, Wilkinson P B, Weller A L, Ogilvy R D, Meldrum P I, and Caunt, S, 2007.
Mineshaft imaging using surface and crosshole 3D electrical resistivity tomography: A
case history from the East Pennine Coalfield, UK. Journal of Applied Geophysics, 62,
324-337.
Chambers, J.E., Gunn, D.A., Wilkinson, P.B., Meldrum, P.I., Haslam, E., Holyoake, S.,
Kirkham, M., Kuras, O., Merritt, A. and Wragg, J., 2014. 4D electrical resistivity
tomography monitoring of soil moisture dynamics in an operational railway
embankment, Near Surface Geophysics, 12, 61-72.
Chivas, A.R. and Nutter, A.H., 1975. Copper Hill porphyry-copper prospect. Knight, S.L. (ed.),
Economic geology of Australia and Papua New Guinea vol. I. Metals: Australasian
Institute of Mining and Meallurgy Monograph No. 5, p. 716-720.
Christensen N.B. and Sorensen K.I., 1994. Integrated use of electromagnetic methods for
hydrogeological investigations. Proceedings of the Symposium on the Application of
Geophysics to Engineering and Environmental Problems, March 1994, Boston,
Massachusetts, 163-176.
Claerbout, J.F. and Muir, F., 1973. Robust modeling with erratic data. Geophysics, 38, 826-
844.
Constable, S.C., Parker, R.L. and Constable, C.G., 1987. Occam’s inversion : A practical
algorithm for generating smooth models from electromagnetic sounding data.
Geophysics, 52, 289-300.
Dahlin, T., 1996. 2D resistivity surveying for environmental and engineering applications. First
Break, 14, 275-284.
Dahlin, T., 2000. Short note on electrode charge-up effects in DC resistivity data acquisition
using multi electrode arrays. Geophysical Prospecting, 48, 181-187.
Dahlin, T. and Bernstone, C., 1997. A roll-along technique for 3D resistivity data acquisition
with multi-electrode arrays, Procs. SAGEEP’97 (Symposium on the Application of
Geophysics to Engineering and Environmental Problems), Reno, Nevada, March 23-26
1997, vol 2, 927-935.
Dahlin, T. and Loke, M.H., 1997. Quasi-3D resistivity imaging-mapping of three dimensional
structures using two dimensional DC resistivity techniques. Proceedings of the 3rd
Meeting of the Environmental and Engineering Geophysical Society. 143-146.
Dahlin,T. and Loke, M.H., 1998. Resolution of 2D Wenner resistivity imaging as assessed by
numerical modelling, Journal of Applied Geophysics, 38, 237-249.
Dahlin, T. and Owen, R., 1998. Geophysical investigations of alluvial aquifers in Zimbabwe.
Proceedings of the IV Meeting of the Environmental and Engineering Geophysical
Society (European Section), Sept. 1998, Barcelona, Spain, 151-154.
Dahlin, T., Bernstone, C. and Loke, M.H., 2002, A 3D resistivity investigation of a
contaminated site at Lernacken in Sweden. Geophysics, 60, 1682-1690.
Dahlin, T. and Zhou, B., 2004, A numerical comparison of 2D resistivity imaging with ten
electrode arrays, Geophysical prospecting, 52, 379-398.
Dahlin, T. and Zhou, B., 2006, Gradient array measurements for multi-channel 2D resistivity
imaging, Near Surface Geophysics, 4, 113-123.
Dahlin, T. and Leroux, V., 2012. Improvement in time-domain induced polarisation data
quality with multi-electrode systems by separating current and potential cables. Near
Surface Geophysics, 10, 545-565.
Dahlin, T. and Loke, M.H., 2015. Negative Apparent Chargeability in Time-domain Induced
Polarisation Data. Journal of Applied Geophysics, in press.
Daily, W. and Owen, E, 1991. Cross-borehole resistivity tomography. Geophysics, 56, 1228-
1235.
Daily, W., Ramirez, A., Newmark, R. and Masica, K., 2004. Low-cost reservoir tomographs
of electrical resistivity. The Leading Edge, May 2007, 272-280.
Daniels F. and Alberty R.A., 1966. Physical Chemistry. John Wiley and Sons, Inc.
deGroot-Hedlin, C. and Constable, S., 1990. Occam's inversion to generate smooth, two-
dimensional models form magnetotelluric data. Geophysics, 55, 1613-1624.
Day-Lewis, F.D., Singha K. , and Binley A., 2005. The application of petrophysical models to
radar and electrical resistivity tomograms: resolution dependent limitations. Journal of
Geophysical Research, 110, B08206, doi: 10.1029/2004JB003569, 17 p.
Dey A. and Morrison H.F. 1979a. Resistivity modelling for arbitrary shaped two-
dimensional structures. Geophysical Prospecting 27, 106-136.
Dey A. and Morrison H.F., 1979b. Resistivity modeling for arbitrarily shaped three-
dimensional shaped structures. Geophysics 44, 753-780.
Edwards L.S., 1977. A modified pseudosection for resistivity and induced-polarization.
Geophysics, 42, 1020-1036.
Ellis, R.G. and Oldenburg, D.W., 1994a, Applied geophysical inversion: Geophysical Journal
International, 116, 5-11.
Ellis, R.G. and Oldenburg, D.W., 1994b. The pole-pole 3-D DC-resistivity inverse problem : a
conjugate gradient approach. Geophys. J. Int., 119, 187-194.
Farquharson, C.G., and D.W. Oldenburg, 1998. Nonlinear inversion using general measures of
data misfit and model structure, Geophysical Journal International, 134, 213-227.
Farquharson, C.G., and D.W. Oldenburg, 2004. A comparison of automatic techniques for
estimating the regularization parameter in non-linear inverse problems. Geophysical
Journal International, 156, 411-425.
Farquharson, C.G., 2008. Constructing piecewise-constant models in multidimensional
minimum-structure inversions. Geophysics, 73, K1-K9.
Fox, R.C., Hohmann, G.W., Killpack,T.J. and Rijo, L., 1980, Topographic effects in resistivity
and induced polarization surveys. Geophysics, 45, 75-93.
Gerard, R. and Tabbagh, A., 1991. A mobile four electrodes array and its application to
electrical survey of planetary grounds at shallow depths. J. Geophys. Res., 96 (B-3),
4117-4123.
Gharibi, M. and Bentley, L.R., 2005. Resolution of 3-D electrical resistivity images from
inversions of 2-D orthogonal lines. Journal of Environmental and Engineering
Geophysics, 10, 339-349.
Golub, G.H. and van Loan, C.F., 1989. Matrix computations. The John Hopkins Un. Press.
Gouet, O., 2007. Projet de classement au titre des articles L.341-1 et suivants du Code de
l’Environnement De L’ÎLE D’OLÉRON. Report, Ministère de l’Écologie et du
Developpement Durable, Département de la Charente-Maritime, Poitiers, 84p.
Griffiths, D.H. and Turnbull, J., 1985. A multi-electrode array for resistivity surveying. First
Break 3 (No. 7), 16-20.
Griffiths D.H., Turnbull J. and Olayinka A.I. 1990, Two-dimensional resistivity mapping with
a computer- controlled array. First Break 8, 121-129.
Griffiths D.H. and Barker R.D.,1993. Two-dimensional resistivity imaging and modelling in
areas of complex geology. Journal of Applied Geophysics, 29, 211-226.
Hallof, P.G., 1990. Reconnaissance and detailed geophysical results, Granite Mountain Area
Pershing County, Nevada. in Fink, J.B., McAlister, E.O., Sternberg, B.K., Wieduwilt,
W.G. and Ward, S.H. (Eds), 1990, Induced polarization : Applications and case
histories : Investigations in Geophysics No. 4, SEG, 325-353.
Harro D. and Kruse S. 2013. Improved imaging of covered karst with the multi-electrode
resistivity implant technique. 13th Sinkhole Conference, NCKRI Symposium 2,
Carlsbad, New Mexico, 7 pp.
Holcombe, J. and Jirack, G., 1984. 3-D terrain corrections in resistivity surveys. Geophysics,
49, 439-452.
Inman, J.R., 1975. Resistivity inversion with ridge regression. Geophysics, 40, 798-817.
Johansson, S. and Dahlin, T., 1996. Seepage monitoring in an earth embankment dam by
repeated resistivity measurements. European Journal of Engineering and Geophysics,
1, 229-247.
Jung, H.K., Min, D.J., Lee, H.S., Oh, S. and Chung, H., 2009. Negative apparent resistivity in
dipole–dipole electrical surveys. Exploration Geophysics, 40, 33–40.
Keller G.V. and Frischknecht F.C.,1966. Electrical methods in geophysical prospecting.
Pergamon Press Inc., Oxford.
Kellett, R.L. and Bauman, P.D., 1999. Electrical resistivity imaging of oil sand horizons : A
modern slant on an old technique. Komex International Ltd.
Kemna, A., Rakers, E. and Binley, A.,1997. Application of Complex Resistivity Tomography
to Field Data from a Kerosene Contaminated Site, In: Proc. 3rd Meeting of the
Environmental and Engineering Geophysics Society, Aarhus, 8-11, September, 1997,
p151-154.
Kenma, A., Binley, A., Ramirez, A. and Daily, W., 2000. Complex resistivity tomography for
environmental applications. Chemical Engineering Journal, 77, 11-18.
Kim, J H., Yi, M. J., Park, S G., Kim, J.G., 2009. 4D inversion of DC resistivity monitoring
data acquired over a dynamically changing earth model. Journal of Applied
Geophysics, 68, 522-532.
Koefoed O.,1979. Geosounding Principles 1 : Resistivity sounding measurements. Elsevier
Science Publishing Company, Amsterdam.
LaBrecque, D.J., Miletto, M., Daily, W., Ramirez, A. and Owen, E., 1996. The effects of noise
on Occam’s inversion of resistivity tomography data. Geophysics, 61, 538-548.
Lee, H., Jung, H.K., Cho, S.H. and Min, D.J. 2014. Negative Apparent Resistivity in 2.5D
Dipole-dipole Electrical Survey for a Very Simple Model. Near Surface Geoscience
2014 - 20th European Meeting of Environmental and Engineering Geophysics.
Li Y. and Oldenburg D.W. 1992. Approximate inverse mappings in DC resistivity problems.
Geophysical Journal International 109, 343-362.
Lines L.R. and Treitel S. 1984. Tutorial : A review of least-squares inversion and its application
to geophysical problems. Geophysical Prospecting, 32, 159-186.
Loke, M.H., 1994. The inversion of two-dimensional resistivity data. Unpubl. PhD thesis, Un.
Of Birmingham.
Loke, M.H., 1999. Time-lapse resistivity imaging inversion. Proceedings of the 5th Meeting of
the Environmental and Engineering Geophysical Society European Section, Em1.
Loke, M.H., 2000. Topographic modelling in resistivity imaging inversion. 62nd EAGE
Conference & Technical Exhibition Extended Abstracts, D-2.
Loke, M.H. and Barker, R.D., 1995. Least-squares deconvolution of apparent resistivity
pseudosections. Geophysics, 60, 1682-1690.
Loke, M.M., 2015. Determination of model reliability in 3-D resistivity and I.P. inversion. 24th
International Geophysical Conference and Exhibition, 15-18 February 2015 - Perth,
Australia.
Loke M.H. and Barker R.D.,1996a. Rapid least-squares inversion of apparent resistivity
pseudosections using a quasi-Newton method. Geophysical Prospecting, 44, 131-152.
Loke M.H. and Barker R.D.,1996b. Practical techniques for 3D resistivity surveys and data
inversion. Geophysical Prospecting, 44, 499-523.
Loke, M.H. and Dahlin, T., 2002. A comparison of the Gauss-Newton and quasi-Newton
methods in resistivity imaging inversion. Journal of Applied Geophysics, 49, 149-162.
Loke, M.H., Acworth, I. and Dahlin, T., 2003. A comparison of smooth and blocky inversion
methods in 2D electrical imaging surveys. Exploration Geophysics, 34, 182-187.
Loke, M.H. and Dahlin, T, 2010. Methods to reduce banding effects in 3-D resistivity
inversion. Procs. 16th European Meeting of Environmental and Engineering
Geophysics, 6 - 8 September 2010, Zurich, Switzerland, A16.
Loke, M.H., Wilkinson, P. and Chambers, J., 2010a. Fast computation of optimized electrode
Detektering av gas i deponier med resistivitet (in Swedish), Rapport SGC 208·1102-
7371, ISRN SGC-R-208-SE, Svenskt Gastekniskt Center, Malmö, 66p.
Roy, A. and Apparao, A, 1971. Depth of investigation in direct current methods. Geophysics,
36, 943-959.
Rucker, D.F., Levin, M.T. and Myers, D.A., 2007. Imaging beneath Hanford's Tank Farms
with electrical resistivity geophysics - An innovative approach. WM '07 Conference,
Feb. 25-March 1, 2007, Tuscon, AZ.
Rucker, D., Loke, M.H., Levitt, M.T. and Noonan, G.E., 2010. Electrical Resistivity
Characterization of an Industrial Site using Long Electrodes. Geophysics, 75, WA95-
WA104.
Rucker, D., Fink, J.B. and Loke, M.H., 2011. Environmental Monitoring of Leaks using Time
Lapsed Long Electrode Electrical Resistivity. Journal of Applied Geophysics, 74, 242-
254.
Rucker, D.F., Crook, N., Glaser, D. and Loke, M.H., 2012a. Pilot-scale field validation of the
long electrode electrical resistivity tomography method. Geophysical Prospecting, 60,
1150-1166.
Rucker, D.F. and Loke, M.H., 2012b. Optimized geoelectrical monitoring with three-
dimensional arrays, Proceedings SEG-AGU Hydrogeophysics Workshop, 8–11 July
2012, Boise, Idaho, USA.
Rucker, D.F. and Noonan, G.E., 2013. Using marine resistivity to map geotechnical properties:
a case study in support of dredging the Panama Canal. Near Surface Geophysics, 11,
625-637.
Rucker, D.F, Crook, N., Winterton, J., McNeill, M., Baldyga, C. A., Noonan, G. and Fink, J.
B., 2014. Real-time electrical monitoring of reagent delivery during a subsurface
amendment experiment. Near Surface Geophysics, 12, 151-163.
Sasaki, Y. 1992. Resolution of resistivity tomography inferred from numerical simulation.
Geophysical Prospecting, 40, 453-464.
Scott, J.B.T., Barker, R.D. Peacock, S., 2000. Combined seismic refraction and electrical
imaging. Procs. 6th Meeting of the European Association for Environmental and
Engineering Geophysics, 3-7 Sept. 1997, Bochum, Germany, EL05.
Seaton, W.J. and Burbey, T.J., 2000. Aquifer characterization in the Blue Ridge physiographic
province using resistivity profiling and borehole geophysics : Geologic analysis.
Journal of Environmental & Engineering Geophysics, 5, no. 3, 45-58.
Schenkel, C.J. and Morrison, H.F., 1994. Electrical resistivity measurement through metal
casing. Geophysics, 59, 1072-1082.
Shima, H., 1992. 2-D and 3-D resistivity imaging reconstruction using crosshole data.
Geophysics, 55, 682-694.
Shima, H., Sakashita, S. and Kobayashi, T., 1996. Developments of non-contact data
acquisition techniques in electrical and electromagnetic explorations. Journal of
Applied Geophysics, 35, 167-173.
Silvester P.P. and Ferrari R.L., 1990. Finite elements for electrical engineers (2nd. ed.).
Cambridge University Press.
Smith, T., Hoversten, M., Gasperikova, E. and Morrison, F., 1999. Sharp boundary inversion
of 2D magnetotelluric data. Geophysical Prospecting, 47, 469-486.
Singer, B.S. and Strack, K.M.,1998. New aspects of through-casing resistivity theory,
Geophysics, 63, 52-63.
Slater, L., Binley, A.M., Zaidman, M.D. and West, L.J., 1997, Investigation of vadose zone
flow mechanisms in unsaturated chalk using cross-borehole ERT. Proceedings of the
EEGS European Section 3rd Meeting, Aarhus, Denmark, 17-20.
Slater, L., Binley, A.M., Daily, W. and Johnson, R., 2000. Cross-hole electrical imaging of a
Wilkinson P.B., Chambers J.E., Meldrum P.I., Ogilvy R.D. and Caunt S., 2006a. Optimization
of array configurations and panel combinations for the detection and imaging of
abandoned mineshafts using 3D cross-hole electrical resistivity tomography. Journal of
Environmental and Engineering Geophysics 11, 213-221.
Wilkinson, P.B., Meldrum, P.I., Chambers, J.C., Kuras, O. and Ogilvy, R.D., 2006b. Improved
strategies for the automatic selection of optimized sets of electrical resistivity
tomography measurement configurations. Geophys. J. Int., 167, 1119-1126.
Wilkinson, P. B., Chambers, J. E., Lelliot, M., Wealthall, G. P. and Ogilvy, R. D., 2008.
Extreme Sensitivity Of Crosshole Electrical Resistivity Tomography Measurements To
Geometric Errors. Geophysical Journal International, 173, 49-62.
Wilkinson, P.B., Loke, M.H., Meldrum, P.I., Chambers, J.E., Kuras, O., Gunn, D.A., Ogilvy,
R.D., 2012, Practical aspects of applied optimized survey design for electrical
resistivity tomography, Geophysical Journal International, 189, 428-440.
Wilkinson, P.B., Uhlemann, S., Chambers, J.C., Meldrum, P. I. and Loke, M.H., 2015.
Development and testing of displacement inversion to track electrode movements on
3D Electrical Resistivity Tomography monitoring grids. Geophysical Journal
International, 200, 1566-1581.
Witherly, K.E. and Vyselaar, J, 1990. A geophysical case history of the Poplar Lake Copper-
Molybdenum deposit, Houston Area, British Columbia. in Fink, J.B., McAlister, E.O.,
Sternberg, B.K., Wieduwilt, W.G. and Ward, S.H. (Eds), 1990, Induced polarization :
Applications and case histories : Investigations in Geophysics No. 4, Soc. Expl.
Geophys.
Wolke, R. and Schwetlick, H., 1988, Iteratively reweighted least squares algorithms,
convergence analysis, and numerical comparisons: SIAM Journal of Scientific and
Statistical Computations, 9, 907-921.
Wynn, J.C., Grosz, A.E. and Carlson, V.L., 1990. Induced-polarization response of some
titanium-bearing placer deposits in Southeastern United States. In Fink, J.B., Sternberg,
B.K., McAlister, E.O., Wieduwilt, W.G. and Ward, S.H. (Eds.), Induced polarization –
Applications and case histories. Investigations in Geophysics No. 4, SEG, 280-303.
Zhao, J.X., Rijo, L. and Ward, S.H., 1986. Effects of geoelogic noise on cross-borehole
electrical surveys. Geophysics, 51, 1978-1991.
Zhdanov, M.S. and Keller, G.V., 1994. The geoelectrical methods in geophysical exploration.
Elseiver, Amsterdam.
Zhou, B. and Greenhalgh, S.A., 1997. A synthetic study on cross-hole resistivity imaging with
different electrode arrays. Exploration Geophysics, 28, 1-5.
Zhou, B. and Greenhalgh, S.A., 2000. Cross-hole resistivity tomography using different
electrode configurations. Geophysical Prospecting, 48, 887-912.
Zhou B. and Dahlin, T., 2003, Properties and effects of measurement errors on 2D resistivity
imaging surveying, Near Surface Geophysics, 1, 105-117.
Zonge, K.L. and Wynn, J.C., 1975. Recent advances and applications in complex resistivity
measurements. Geophysics, 40, 851-864.
Zhu, T. and Feng, R., 2011. Resistivity tomography with a vertical line current source and its
applications to the evaluation of residual oil saturation. Journal of Applied Geophysics,
73, 155–163.
In this equation αs is the relative weight given to the constraint such that the inversion model q
is close to a reference model q0. Note it is basically the damping factor used in the Marquardt-
Levenburg method (Line and Treitel, 1984). Normally a small value (such as 0.05) is used for
the relative damping factor weight αs. In this case the smoothness constraint in the FR roughness
matrix is the dominant constraint.
pole-dipole array with a constant ‘a’ spacing as well as the generally poorer resolution of the
I.P. data. The other reason is that the use of the smoothness constraint tends to produce a model
with slowly varying model values. The I.P. data does not have much information in the lower
part of the model. The high I.P. values towards the bottom are artifacts due to the use of the
smoothness constraint. Increasing the relative damping factor weight αs to 0.20 cause slight
charges in the I.P. model (Figure 156b). Using a higher value of 0.50 starts to close off the
bottom of the I.P. anomaly (Figure 156c). The constraint associated with αs forces model to be
‘close’ to a homogenous background model with a value of close to 0 mV/V. Increasing αs
further to 1.0 further improves the shape of I.P. anomaly that is closer to the true structure
(Figure 156d). In this case, equal weights are given to the model smoothness and deviations
from the constant background model. The model I.P. values below the block are reduced as
there is not much data to support higher values there and the effect of the smoothness constraint
is reduced. The same value for damping factor λ was used in all the inversions. It should be
possible to improve the results, particularly for the larger values of αs, by using an adaptive
method such as the L-curve (Farquharson and Oldenburg 2004) to automatically set this
damping factor.
Figure 156. Models obtained with different relative damping weights αs.
Figure 157. Models for the Burra data set with a relative damping weight of 0.5 for αs.
Figure 158. I.P. vertical sections along the y-direction (at x=1650 to 1700 m) Burra data set
(a) without (αs=0.0) and (b) with (αs=0.5) a reference model constraint.
Figure 159. A long electrode in a (a) homogeneous medium, (b) two-layer medium with a low
resistivity lower layer and (c) partly in air.
C.1 The distributed receivers I.P. systems and negative apparent resistivity values.
In most surveys, the measured apparent resistivity values (geometric factor multiplied
by the measured resistance) are always positive. Negative apparent resistivity values are
usually caused by noise in the measurements. However, in some situations, negative apparent
resistivity values have been observed which are not due to noise (Lee at el., 2014). It has
become more common with the use of the offset type of arrays (White et al., 2001) and non-
conventional arrangements (Figure 160) used in large scale 3-D I.P. surveys, particularly with
near-surface large resistivity contrasts and topography. In I.P. surveys, the circuit for the
current electrodes at the transmitter is always separate from the receiver dipoles. In the older
I.P. systems, it is basically impossible to determine whether the measured potential has the
same sign as the transmitted current. Thus, a crude approach has been to assume the apparent
resistivity value is always positive. This causes a problem in modeling the data when the finite-
difference or finite-element program correctly predicts that some measurements should be
negative. It is impossible to determine whether the input apparent resistivity is wrong, or it is
due to numerical errors in the forward modeling method! Fortunately, the newer I.P. systems
(such as the Iris Instruments Full-Waver system) now has synchronization features so that it
possible to determine the correct sign of the apparent resistivity values. It has been found that
5 to 10% of the measurements might have negative apparent resistivity values, particularly with
the distributed type of system (Figure 160b), and it is important to correctly identify the
negative apparent resistivity values for a more accurate resistivity inverse model. Note the sign
of the apparent I.P. values can always be correctly determined as it is measured relative to the
primary resistivity potential.