- Email: [email protected]
The effect of dike geometry on different resistivity configurations
Journal of Applied Geophysics 57 (2005) 278 – 292 www.elsevier.com/locate/jappgeo
The effect of dike geometry on different resistivity configurations Thomas Henniga,T, Andreas Wellera, Tran Canhb a
Institut fu¨r Geophysik, Technische Universita¨t Clausthal, Arnold-Sommerfeld-Str. 1, D-38678 Clausthal-Zellerfeld, Germany b Institute of Geological Sciences, Vietnamese Academy of Science and Technology, Hanoi, Vietnam Received 19 March 2004; accepted 1 March 2005
Abstract Geoelectrical profiling with multi-electrode systems has become an important tool for monitoring dike embankments bordering rivers. Profiles running perpendicular to the dike axis are affected by the dike topography, with the amplitude of this effect dependent on the surface geometry and the choice of the electrode configuration. Investigations using seven different electrode configurations have shown that some configurations are less sensitive to the topography than others. The topography correction method (TCM) is an important tool for processing data from measurements at river dikes. This method is generally recommended for flank angles steeper than 108. The topography effect is calculated by two-dimensional finite element modelling. The resulting synthetic data of a homogeneous dike body are used to apply a topographic correction for each measurement. The topographic effect and correction procedure is demonstrated for synthetic dike data and for a data set from a river dike in Thai Binh province (Vietnam). The topography can be ignored for flank angles less than 258 if an averaged Half-Wenner electrode configuration is used. This configuration has proved to be less affected by undulated topography and the focusing effect of averaging the two data sets provides reliable structural information without the need for time-consuming data inversion. D 2005 Elsevier B.V. All rights reserved. Keywords: Electrical resistivity; Finite element method; Topography correction; Dike monitoring
1. Introduction Multi-electrode arrays are widely used in monitoring dike embankments that border rivers to prevent flooding (Ullrich and Meyer, 2003). Cables
T Corresponding author. Fax: +49 5323 722320. E-mail address: [email protected] (T. Hennig). 0926-9851/$ - see front matter D 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.jappgeo.2005.03.001
connect a large number of equally spaced electrodes with a control unit. The active electrodes are selected by the control unit according to a predefined measurement sequence. A set of different electrode configurations which are commonly used in multi-electrode surveys is shown in Fig. 1 with A and B being the current electrodes, M and N the potential electrodes, a the fixed electrode distance and n an integer factor. Each configuration has its
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
279
(a) Dipole - Dipole B
A
M
a
N
na
a
(b) Wenner-α A
M
N
na
na
B na
(c) Wenner-β B
A
M
na
na
N na
(d) averaged Wenner-α/β A
M
N na
na B
A
B na
M
N
(e) Half-Wenner forward A
M na
N
B
oo
A
B
oo
B
oo
na
(f) Half-Wenner backward N
M na
na
(g) averaged Half-Wenner A
M
N
na N
na M
A
(h) Schlumberger A
M na
N a
B na
Fig. 1. Commonly used electrode configurations for a multi-electrode survey. (a) Dipole–dipole. (b) Wenner-a. (c) Wenner-h. (d) Averaged Wenner a/h. (e) Half-Wenner forward. (f) Half-Wenner backward. (g) Averaged Half-Wenner. (h) Schlumberger.
advantages and drawbacks (Ward, 1990). The Wenner-a (Fig. 1b) configuration is suitable to resolve resistivity changes with depth. Gradient arrays like dipole–dipole configuration (Fig. 1a) provide a better lateral resolution but they become
severely affected by noise as the electrode spacing increases. Earlier investigations have shown that a special type of three electrode configuration which is called Half-Wenner configuration (Fig. 1e,f) provides some
280
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
advantages for resistivity profiling with multi-electrode systems (Peschel, 1967): (i) The number of measurements with a fixed number of electrodes is larger than for Wenner-a and dipole–dipole configuration. Since more levels can be measured an increased depth of penetration can be reached. (ii) For the infinite electrode B a good location for current injection can be found a large way outside the profile. (iii) Two Half-Wenner readings can be combined to a Wenner-a reading. (iv) The simple combination of a Half-Wenner forward configuration and a Half-Wenner backward configuration with a fixed position of electrode M has proved to act as a focussing tool (Kampke et al., 1998). The average of both readings, which is called the averaged HalfWenner (Fig. 1g), results in a pseudo-section that shows the main features of the subsurface resistivity distribution. It will be shown in this paper, that the averaged Half-Wenner configuration is also less sensitive than other commonly used configurations (e.g. Wenner-a, dipole–dipole) to a rough topography. Resistivity profiling is normally applied along the dike axis. Profiles are measured either at the top of the dike, the flanks, or even at the foot of the dike to explore the foundation. If some resistivity anomalies are found it is often difficult to predict the exact location of the anomalous structure. A two-dimensional interpretation assumes that the anomalies are caused by structures directly beneath the profile. But it is known that structures at both sides of the profile also affect the measurements. In order to find out the exact location of the zone of anomalous resistivity it is advisable to measure a perpendicular profile, which provides a cross-section of the dike body. Such a profile will usually be severely affected by the dike geometry. Different authors have investigated the influence of topography on resistivity readings (e. g. Fox et al., 1980; Tsourlos et al., 1999). Fox et al. (1980) have shown that a dipole–dipole survey crossing a ridge of homogeneous resistivity causes an apparent resistivity maximum under the top and minima at the two feet of the ridge. The maximum exceeds the intrinsic
resistivity by the factor 2.6. The minimum drops to 50% of the true resistivity. Tsourlos et al. (1999) used a topography similar to the dike geometry. The resulting resistivity pseudo-section for a dipole–dipole array shows a strong resistivity maximum under the dike body and adjacent minima near the feet of the dike. The values of the maximum and the minima are comparable to the model of Fox et al. (1980) that uses the same slopes of 308 for both flanks. Tsourlos et al. (1999) stated that slopes larger than 108 cause misleading artificial errors that have to be removed. In this paper, the effect of topography for different electrode configurations is demonstrated for both a model of synthetic data and field data from a Vietnamese dike. A topography correction method is applied (Fox et al., 1980) that largely eliminates the influence of the dike geometry from the measured data.
2. Modelling with finite element method (FEM) Forward modelling of geoelectrical measurements can be done by different numerical methods: boundary element method (e.g. Dieter et al., 1969; Okabe, 1981; Schulz, 1985), finite difference modelling (e.g. Mufti, 1976; Dey and Morrison, 1979a,b; Weller et al., 1996) or finite element method (FEM). Coggon (1971) was one of the first authors to use FEM for two-dimensional problems. An overview can be found in Hohmann (1988). The FEM provides some advantages for modelling geoelectrical data for an environment with a rough topography as the triangular elements are more flexible for more complicated structures compared with the rectangular elements of finite difference modelling. The Neumann-type boundary condition at the earth’s surface is automatically fulfilled by the formulation of the variation problem. Additionally, the discretisation errors at the surface are smaller than the errors using finite difference modelling due to the correct implementation of the boundary plain between airspace and subsurface for undulated topography. The actual problem of geoelectrical modelling is three-dimensional (3-D). The electrical potential U(x,y,z) in a non-uniform isotropic medium is governed by the following differential equation j½rð x; y; zÞU ð x; y; zÞ þ jjs ð x; y; zÞ ¼ 0;
ð1Þ
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
where the function r(x,y,z) describes the conductivity distribution and js(x,y,z) the source distribution of current density. A Neumann-type boundary condition at the surface and asymptotic conditions at the other model boundaries complete the problem. The differential Eq. (1) is valid for steady current flow and is also acceptable for alternating current at low frequencies if electromagnetic effects can be ignored. Dike structures can be represented by 2-D models. If the y-axis of the model is chosen to be parallel to the strike direction, the conductivity distribution becomes independent of y and can be written as r(x,y,z) = r(x,z). Even a 2-D resistivity structure causes a 3-D potential distribution. A pure 2-D modelling requires line sources in the strike direction (Mufti, 1976). However, most of the current electrodes used in geoelectrical prospecting can be assumed to be point sources, and the problem remains 3-D. The numerical expense in solving such problems can be reduced by using a Fourier-cosine transformation of potentials in strike direction y (Dey and Morrison, 1979a): Z l ˜ U ð x; k; zÞ ¼ 2 U ð x; y; zÞcosðkyÞdy; ð2Þ 0
where k denotes the wavenumber. The current source is assumed to be located in the plane y = 0. The application of the Fourier-cosine transformation to the 3-D Eq. (1) results in a 2-D partial differential equation: B B B B r U˜ þ r U˜ Bx Bx Bz Bz ð3Þ 2 ˜ rk U þ Idð x xs Þdð z zs Þ ¼ 0; where I is the current which is injected at position (x s,0,z s) and d is the Dirac delta function. Eq. (3) has to be solved for a set of discrete wavenumbers k. The associated boundary condition for the upper boundary is a Neumann-type condition with r air = 0. A mixed boundary condition as proposed by Dey and Morrison (1979b) is applied at the left, right and lower borderline: rndjU bU ¼ 0;
ð4Þ
with n being the unit vector perpendicular to the boundary and b a factor that defines the asymptotic behaviour of the potential field in a homogeneous
281
medium. Compared with the often used Dirichlet condition, it allows a higher accuracy with a smaller modelling area (Weller et al., 1996). Consequently, fewer elements and less calculating time are required. The discretisation of the subsurface is shown in Fig. 2. This model represents an idealized dike with flank angles of 308. The highest density of nodes is chosen in the inner area around the electrodes. The modelling area is extended to the left and right and to the bottom to satisfy the mixed boundary conditions. The inter-nodal distance increases continuously outside the inner area. Each quadrangular element is cut into four triangles. This implements an additional node in the centre of the element and consequently a higher level of discretisation of the model area. Using the static condensation procedure (Schwarz, 1991), the additional nodes are eliminated but the higher level of discretisation is kept. Thus, a better discretisation can be attained without increasing the dimension of the system of equations, while the calculation time is reduced. ˜ (x,k,z) is assumed to The transformed potential U vary linearly over each of the triangular elements. Combining the resulting equations for all elements yields a symmetric system of equations with a dimension equal to the number of nodes. A characteristic property of FEM is that the resulting matrix is symmetric, sparse, banded and diagonally dominant. The structure depends on the order in which the nodes are numbered. The bandwidth of the resulting matrix, which is equal to the maximal index difference of nodes in all elements, should be minimized. The properties of the matrix are used to solve the linear system of equations with the Cholesky algorithm (Schwarz, 1991) for all wavenumbers. The resulting ˜ (x,k,z) in the wavenumber domain are potentials U used to reconstruct the potential U(x,y,z) by an inverse Fourier-cosine transform. If the potential electrodes are located in the plane y = 0 the inverse Fourier-cosine transformation is performed by a simple integration Z 1 l ˜ U ð x; y; zÞ ¼ U ð x; k; zÞdk; ð5Þ p 0 which is done numerically by a combination of Gaussian quadrature and Laguerre integration as proposed by LaBreque et al. (1996). The choice of the wavenumbers affects the accuracy and the computing time of the algorithm. In our algorithm, nine wave
282
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292 4
2
electrode position finite element grid node
z in m
0
-2
-4
-6
-8 -4
-2
0
2
4
x in m
6
8
10
12
14
Fig. 2. Ridge model for synthetic modelling after Fox et al. (1980) with the discretisation used for FE modelling.
numbers are used in order to keep the error of the inverse transformation less than 0.5% (Kemna, 2000). For each electrode configuration the potential difference between the electrodes M and N and the geometry factor K are used to calculate the apparent resistivity for each measurement: qa ¼
UM UN K: I
ð6Þ
The forward problem is solved for each current source separately. The apparent resistivities for all relevant electrode configurations can be computed easily by a superposition of several pole–pole configurations. This method reduces the number of forward problems to the number of current electrodes. The number of measured configurations is generally much higher.
3. Topographic correction method (TCM) Topographic effects occur if the surface is not flat. The intensity of the effect depends first of all
on the roughness of the surface and additionally on the electrode configuration used. Assuming a flat surface with a homogeneous subsurface, a current injected at an electrode at the surface causes a potential field around the electrode where the isopotential planes are represented by concentric halfspheres in the subsurface. According to the Neumann-type boundary condition no current can enter into the air (r air = 0). Consequently, the potential lines are perpendicular to the surface. Considering an uneven topography, the potential planes have to act in the same way, they remain perpendicular to the undulating surface. Consequently, the subsurface potential distribution is disturbed. If this effect is ignored anomalies in the subsurface would not be correctly interpreted. In order to determine the topographic effect of a certain configuration the potential distribution for any current electrode position is calculated, using a constant resistivity q mod, for all elements of the model. The apparent resistivity q cal results from the potential difference between the positions of electrodes M and N using the flat surface geometric factor
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
K that ignores the small changes in electrode distances caused by the topography. In the case of a flat surface for all configurations the apparent resistivity is expected to be close to the constant resistivity q mod attributed to all elements. Slight differences, that are called discretisation errors, are caused due to the finite accuracy of the finite element algorithm. Assuming an undulating surface, the resulting apparent resistivity q cal may differ considerably from the intrinsic value q mod. In order to correct the resulting pseudo-sections for the topographic effect correction factors CF (Fox et al., 1980) are determined by dividing the constant resistivity value q mod, that is attributed to all elements, by the calculated apparent resistivity q cal: CF ¼
qmod : qcal
ð7Þ
The CF that are calculated for a given surface topography can then be applied to the measured apparent resistivity data. Only if surface topography alone causes a resistivity anomaly the TCM yields exactly the data that would have been measured on a flat surface. Considering both inhomogeneous resistivity structures and an undulating surface, the topographic correction by the CF assumes that topography response and the signal resulting from the resistivity target are linearly superimposed. Since the resistivity forward problem is nonlinear this assumption can only be an approximation. But comprehensive tests with simple resistivity structures and surface geometry have shown that the TCM yields satisfactory results (Fox et al., 1980; Tsourlos et al., 1999).
4. Synthetic modelling experiments 4.1. Comparison with the ridge model The developed algorithm is tested in comparison with the results of Fox et al. (1980). Therefore, an almost identical ridge model is used to determine the apparent resistivity for a homogeneous subsurface. The grid used is shown in Fig. 2. It consists of 405 nodes and 370 elements. The calculated results match the results of Fox et al. (1980) for the dipole–dipole electrode configuration. Some deviations are found in
283
the shallow subsurface. This effect probably occurs because of the different discretisation of the subsurface air interface in both algorithms. Fox et al. (1980) use finite element technique but they fill up the modelling area to a horizontal line above the surface with elements of high resistivity values. Though Fox et al. (1980) do not exactly describe the dimension of the elements used, the average difference between the published values and our results reaches eight per cent. The largest differences (up to 20%) appear for the configurations with minimal spacing. Further tests with the dike model from Tsourlos et al. (1999) show a slightly better agreement with our modelling results. The average difference reaches 7% (Hennig, 2003). 4.2. Dike model with termite nest River dikes in the northern provinces of Vietnam have a total length of about 5000 km. The social and economic development of this region largely depends on the integrity of the dike system to prevent flooding. The safety of river dikes is affected by a diversity of problems. Sandy formations in the dike foundations cause seepage effects. Various termite species dig cavities for their nests in the dike body and decrease the stability of the dike (Tuyen et al., 2000). The efficiency of geophysical methods for dike monitoring has to be improved to detect such defects in the dikes. To study the problem of termite nests, an asymmetric dike model with a resistive cavity has been investigated. The dike model which is shown in Fig. 3 consists of a steeper flank with a slope of 308 on the left side and the other flank descends with an angle of 158. The cavity with a cross-section area of 4 m2 is placed 1.5 m below the top of the dike. It represents a typical termite nest with a resistivity value of 1000 V m (black box in Fig. 3). A resistivity of 25 V m is attributed to the surrounding silty dike material. The grid for FEM consists of 3941 nodes and 3756 elements. Two apparent resistivity data sets of a Wenner-a configuration have been generated by FEM modelling. The first data set corresponds to the homogeneous dike cross-section without a termite nest. The resulting pseudo-section which is displayed in Fig. 4a is shifted vertically according to the dike surface geometry. The depth of penetration is assumed to be
284
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
15
z in m
10
termite nest 1000 Ωm
5 dikebody 25 Ωm
0 -5 15
20
25
30
35
40
45
50
55
60
65
70
x in m Fig. 3. Dike model for synthetic modelling of a river dike with an embedded resistive structure (1000 V m) in the dike body (25 V m).
(a) z in m
10 5 0
15
20
25
30
35
40
45
50
55
60
65 31
x in m
30 29
(b)
28 10
z in m
27 26
5
24 0
23 22 15
20
25
30
35
40
45
50
55
60
65
21
x in m
20 19
(c) ρ in
z in m
10
m
5 0
15
20
25
30
35
40
45
50
55
60
65
x in m Fig. 4. Resistivity pseudo-sections of synthetic Wenner-a data of dike model in Fig. 3. (a) Model without termite nest, pseudo-section without topographic correction. (b) Model with termite nest, pseudo-section without topographic correction. (c) Model with termite nest, pseudo-section after application of topographic correction.
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
na/2 for all pseudo-sections. This is a good approximation for Wenner configurations (Barker, 1989). Though no nest was modelled, the pseudo-section shows a resistive structure at shallow depth beneath the top of the dike. But this anomaly and the other resistive and conductive features result only from the dike topography. The second data set was generated including the nest in the dike body. The resulting pseudo-section (Fig. 4b) shows a more accentuated resistive anomaly at the location of the nest. But it becomes obvious from the comparison of both pseudo-sections that a reliable conclusion whether a nest is located in the dike can not be drawn. The application of TCM to the second data set according to the previously described procedure yields the corrected pseudo-section which is shown in Fig. 4c. Though the lateral resolution of Wenner-a configuration is limited, an extended weak anomaly with resistivity values up to 28 V m becomes visible beneath the top of the dike. A similar anomaly would have been generated from a cavity below a flat surface. Most other anomalous structures which were related to the dike topography disappear after application of TCM. It should be noted that the application of TCM to the synthetic pseudo-section in Fig. 4a results in a homogeneous section with constant resistivity of 25 V m. Another data set that consists of both forward and backward Half-Wenner data has been generated by FEM modelling. A 2-D inversion of this data set was performed with the program DC2DSIRT (Kampke, 1996) which is based on a finite difference forward procedure (Weller et al., 1996) and a simultaneous iterative reconstruction technique. This robust inversion technique, which uses a back-projection as starting model, yields smooth models accentuating the main features of subsurface resistivity distribution after few iterations. The software which enables the use of arbitrarily chosen configurations has been successfully applied for hydrogeological (Olayinka and Weller, 1997) and archaeological (Kampke, 1999; Schleifer et al., 2002) projects. Several tests have proved that the results are comparable with those of commercial programs (e.g. Loke and Barker, 1995). Fig. 5a represents the inversion result after ten iterations. A resistive anomaly can be found in the deeper subsurface with a maximum value of approximately 50 V m. Almost no indication can be detected
285
in the shallow subsurface region of the termite nest. The topographic effect hides the target and causes a misleading interpretation. The termite nest would not have been detected. Applying TCM to the modeled data set and inverting the corrected data set yields the cross-section shown in Fig. 5b. A resistive anomaly is situated exactly at the predetermined position. The deeper structures, which were caused by the dike topography, disappear completely. Resistivity values close to 25 V m are calculated for the dike material. The resistive anomaly that represents the termite nest reaches about 30 V m. Though the maximal value does not correspond to the true value of 1000 V m, the indicated structure shows the termite nest in the correct position. These examples show the influence of topography on geoelectrical data, demonstrating that the interpretation of uncorrected data sets can be wrong, especially for small anomalies like termite nests. The data quality in the theoretical model is ideal whereas field measurements are often affected by noise. Consequently, TCM becomes more important for field data. 4.3. Sensitivity to topography of different electrode configurations The topographic surveying results of a dike at the right side of Tra Li river in Thai Binh province have been selected to demonstrate the sensitivity of different electrode configurations to a typical dike geometry. The dike cross-section and the FEM grid are shown in Fig. 6. A homogeneous resistivity distribution is assumed in the whole dike body and the subsurface. The height of the dike is about 10 m and the length of the profile crossing the dike is about 35 m. The flank angles reach up to 238. The correction factors (CF) are calculated for seven different configurations: ! ! ! ! ! ! !
dipole–dipole (Fig. 1a), Wenner-a (Fig. 1b), Wenner-h (Fig. 1c), averaged Wenner-a/h (Fig. 1d), Half-Wenner forward (Fig. 1e), Half-Wenner backward (Fig. 1f), and averaged Half-Wenner (Fig. 1g).
286
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
(a) z in m
10 5 0
15
20
25
30
35
40
45
50
55
60
65
x in m ρ in
m 18
22
26
30
34
38
42
46
50
54
(b) z in m
10 5 0
15
20
25
30
35
40
45
50
55
60
65
x in m ρ in
m 24
25
26
27
28
29
30
31
Fig. 5. Inversion results of synthetic Half-Wenner forward and backward data of dike model in Fig. 3. (a) Without topographic correction. (b) With topographic correction.
A summary of all CF and the total number of data are compiled in Table 1. The resulting pseudo-sections of CF are embedded in the dike geometry plots of Figs. 7 and 8. The influence of topography on the
resistivity readings can be seen in all configurations. However, the intensity and spreading of CF-values varies over a wide range for the considered configurations. The commonly used dipole–dipole config-
electrode position
z in m
5
finite element grid node
0
-5
-10 -25
-20
-15
-10
-5
0
x in m
5
10
15
20
Fig. 6. Topographic profile of river dike in Thai Binh province with discretisation.
25
30
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
287
Table 1 Values of correction factors (CF) for different electrode configurations for the dike geometry shown in Fig. 6 Dipole–dipole Wenner-a Wenner-h Half-Wenner forward Half-Wenner backward Averaged Half-Wenner Wenner-a/h averaged
Number of data
Minima of CF
Maxima of CF
Mean of CF
S.D. of CF
904 664 664 800 800 1600 1328
0.69 0.84 0.65 0.68 0.72 0.90 0.93
1.43 1.30 1.27 1.24 1.41 1.10 1.10
1.01 1.01 0.98 0.98 1.01 1.00 1.00
0.13 0.08 0.14 0.13 0.16 0.05 0.04
uration (Fig. 7a) is highly affected by topography. The CF values range from 0.69 to 1.43. According to Eq. (7) a CF of 1.43 corresponds to 43% deviation of the
apparent resistivity q cal from the true resistivity q mod. Beside the significant topography effect dipole–dipole data are generally more disturbed by noise, especially
(a) z in m
5
0
-10
-5
0
5
10
15
-10
-5
0
5
10
15
-10
-5
0
5
10
15
-10
-5
0
5
10
15
(b) z in m
5
0
(c) z in m
5
0
(d) z in m
5
0
correction factor
0.6
0.8
x in m
0.95
1.1
1.3
1.5
1.7
Fig. 7. Correction factors for different electrode configurations. (a) Dipole–dipole. (b) Wenner-a. (c) Wenner-h. (d) Averaged Wenner-a/h.
288
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
(a) z in m
5
0
-10
-5
0
5
10
15
-10
-5
0
5
10
15
-10
-5
0
5
10
15
(b) z in m
5
0
(c) z in m
5
0
correction factor 0.6
0.8
x in m
0.95
1.1
1.3
1.5
1.7
Fig. 8. Correction factors for different electrode configurations. (a) Half-Wenner forward. (b) Half-Wenner backward. (c) Averaged HalfWenner.
for larger dipole distances n. Therefore, the depth of penetration is normally restricted. The spreading of CF-values is slightly lower for Wenner-h (Fig. 7c), and the Wenner-a configuration yields a further reduction, where the CF-values range from 0.84 to 1.30 (Fig. 7b). Similar CF variations were obtained from the HalfWenner forward configuration that is shown in Fig. 8a. The Half-Wenner backward CF vary between 0.72 and 1.41 (Fig. 8b). Both Half-Wenner forward and backward configurations are also sensitive to topography. On closer examination of both diagrams it becomes obvious that the minima of the forward mode coincide with the maxima of the backward mode and vice versa. Consequently, the average of both images (Fig. 8c) results in considerably smaller CF varying between 0.91 and 1.10. Thus, the sensitivity to topography decreases for the averaged configuration compared to the original forward and backward modes.
The averaged Wenner-a and Wenner-h (Fig. 7d) as well as the averaged Half-Wenner configuration (Fig. 8c) are less sensitive to rough topography. The CFvalues vary from 0.90 to 1.10 for averaged Wenner-a and Wenner-h. As shown in Table 1, not only the interval but also the standard deviation (S.D.) of the CF-values becomes significantly reduced for the averaged configurations. The standard deviation can be regarded as topography induced noise that is added to the measured signal. Since the suppression of noise is a key criterion of data acquisition, configurations with low standard deviation in CF-values should be preferred. As mentioned above, the number of measurements for the same number of electrodes is higher for HalfWenner configurations (1600) than for Wenner-a/h (1328) in this example. The information density of the subsurface becomes more than 15% higher. Regarding this result and the above-mentioned advantages, it is recommended to use the averaged
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
289
device (GeoSys Leipzig) in Half-Wenner forward and backward electrode configurations. The profile length is 33 m. A constant electrode spacing of 0.5 m was used along the profile. The maximum level of investigation is n = 16. The uncorrected Half-Wenner forward, Half-Wenner backward, Wenner-a and averaged Half-Wenner pseudo-sections are presented in Fig. 9. The Wenner-a data were generated by averaging two Half-Wenner readings with the same electrode spacing na and a common reference point in the centre between the electrodes M and N. It can be seen that the asymmetric configurations, like HalfWenner forward and backward, result in asymmetric
Half-Wenner configuration for geoelectrical data acquisition. The vertical and lateral resolution of this averaged configuration satisfies the requirements of dike investigations.
5. Field data The effect of topographic correction is demonstrated in Figs. 9 and 10 for a field data set from Tra Li river in Thai Binh province (Vietnam). The dike topography is the same as in the example shown in Fig. 6. The data was acquired with a GSM 150
(a) z in m
5
0 -10
-5
0
5
10
15
20
-10
-5
0
5
10
15
20
-10
-5
0
5
10
15
20
-10
-5
0
5
10
15
20
(b) z in m
5
0
(c) z in m
5
0
(d) z in m
5
0
x in m a
in
m 12
18
26
37
55
80
Fig. 9. Pseudo-sections of field data from Tra Li river dike. (a) Half-Wenner forward. (b) Half-Wenner backward. (c) Wenner-a. (d) Averaged Half-Wenner.
290
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
(a) z in m
5
0
-10
-5
0
-10
-5
0
5
10
15
20
5
10
15
20
z in m
(b)5
0
x in m in
m 12
18
26
37
55
80
Fig. 10. Inversion result of field data from Tra Li river dike. (a) Without TCM. (b) With TCM.
pseudo-sections. Some maxima of the forward measurement correspond with minima of the backward measurement and vice versa. As shown in the previous chapter, this topographic effect is widely reduced by an averaging overlay of Half-Wenner forward and backward pseudo-section. Compared with the Wenner-a pseudo-section, which is more distorted by topography, the averaged Half-Wenner image is smoother. The measured data set, consisting of 1148 HalfWenner forward and backward readings, was processed with the 2-D inversion program DC2DSIRT. In a first run the original data were inverted for a flat surface, whereas for the second run the topography correction was applied to the data. The resulting resistivity models are shown in Fig. 10a and b. Since only slight differences between the two images are visible it becomes obvious that a topography correction can be neglected if both Half-Wenner forward and backward data are available. The comparison with the averaged pseudo-section in Fig. 9d verifies that the inversion does not change the image significantly. The high imaging quality of the averaged Half-Wenner pseudo-section results from the focusing effect of this configuration (Kampke, 1999). All relevant features become visible in both the pseudo-section and the inverted model. A resistive anomaly is found in the shallow subsurface of the land side flank of the dike at
x = 2 m. As an excavation shortly after the survey has proved, this anomaly is caused by an extended termite nest. Since no excavations were allowed at the river side dike flank, the origin of the other resistive anomaly at x = 7 m could not be verified.
6. Conclusions Geoelectrical multi-electrode measurements are a suitable tool for investigating river dikes although topographic effects often complicate the interpretation of the measured data. To overcome the problems caused by a rough topography either topographic correction methods (TCM) or less sensitive electrode configurations should be used. TCM is based on a finite element forward modelling that considers the real dike geometry. The dike body is assumed to be homogeneous. The resulting section of correction factors (CF) represents the topographic effect at the measured data. The TCM uses CF to remove topographic effects from measured data. The interval of CF-values and their standard deviation enable a quantitative comparison of the sensitivity to topography of different electrode configurations. It has been shown that an appropriate combination of different electrode configurations reduces the
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
sensitivity of the generated pseudo-sections to the topography. This effect was demonstrated with the averaged Half-Wenner configuration and the average of the Wenner-a and Wenner-h configurations. The investigations have shown that the topographic effect becomes larger as the flank angle increases. The topographic correction can be neglected for both averaged configurations if the flank angles are less than 258. Since this condition is fulfilled at most dikes there is a clear advantage of using these configurations for dike monitoring. Other configurations that are more affected by topography like dipole–dipole and Wenner-a require a topographic correction even at lower flank angles. The use of configurations that are less affected by topography makes the interpretation also less sensitive to uncertainties of the geometric surveying data. The larger the topographic effects of the chosen configuration the more exact the surface topography has to be surveyed. Thus, the expense for geometric surveying can be reduced if appropriate configurations are applied. This conclusion remains valid even if modern inversion procedures which include an automatic topographic correction are used. Another advantage of the averaged configurations is the resemblance between the pseudo-section and the section obtained by inversion that results from the focusing effect of these combined configurations. Without any time-consuming inversion a reliable image of the subsurface structures can be generated in the field. Termite nests and other defects in the dike body can be easily detected in the pseudo-section plot. The operator can quickly provide the relevant information that would allow more detailed lines to be measured right away. Considering the low sensitivity to topography, the high image quality and the good data coverage, a survey with a combination of Half-Wenner forward and backward configurations is highly recommended for dike monitoring.
Acknowledgments The presented investigations were made in the joint scientific project bDevelopment of a geoelectrical technology to find termite nests in river dikesQ
291
between the Institute of Geological Sciences at the Vietnamese Academy of Sciences and Technology and the Institute of Geophysics at Technical University Clausthal (Germany) that was sponsored by the Volkswagen-Stiftung (2001–2004). We thank for the assistance provided by the authorities of Thai Binh province which enabled our field measurements at a dike section at Tra Li river.
References Barker, R.D., 1989. Depth of investigation of collinear symmetrical four-electrode arrays. Geophysics 54, 1031 – 1037. Coggon, J.H., 1971. Electromagnetic and electrical modeling by the finite element method. Geophysics 36, 132 – 155. Dey, A., Morrison, H.F., 1979a. Resistivity modelling for arbitrarily shaped two-dimensional structures. Geophys. Prospect. 27, 106 – 136. Dey, A., Morrison, H.F., 1979b. Resistivity modeling for arbitrarily shaped three-dimensional structures. Geophysics 36, 753 – 780. Dieter, K., Paterson, N., Grant, F., 1969. IP and resistivity type curves for three-dimensional bodies. Geophysics 34, 615 – 632. Fox, R.C., Hohmann, G.W., Killpack, T.J., Rijo, L., 1980. Topographic effects in resistivity and induced polarization surveys. Geophysics 45, 75 – 93. Hennig, T., 2003. Modellierung von Polarisationseffekten mit finiten Elementen. Diplomarbeit. Institut fqr Geophysik, Clausthal-Zellerfeld. Hohmann, G.W., 1988. Numerical modeling for electromagnetic methods in geophysics. In: Nabighian, M.N. (Ed.), Electromagnetic Methods in Applied Geophysics, vol. 1. Invest. Geophys. Soc. Exp. Geophys, Tulsa, OK, pp. 313 – 363. Kampke, A., 1996. Modellierung and Inversion von Daten der induzierten Polarisation fqr zweidimensionale Verteilungen der elektrischen Leitf7higkeit. Diplomarbeit. Institut fqr Geophysik und Meteorologie, Braunschweig. Kampke, A., 1999. Focused imaging of electrical resistivity data in archaeological prospecting. J. Appl. Geophys. 41, 215 – 227. Kampke, A., Weller, A., Peschel, G., 1998. Focussing effect on an averaged three-electrode configuration. Proceedings of the IVth Meeting of EEGS-European Section, Barcelona, pp. 865 – 868. Kemna, A., 2000. Tomographic inversion of complex resistivity— theory and application. Bochum, Ruhr-Universit7t, Dissertation. LaBreque, D.J., Miletto, M., Ramirez, W., Owen, E., 1996. The effects of noise on Occam’s inversion of resistivity tomography data. Geophysics 61, 538 – 548. Loke, M.H., Barker, R.D., 1995. Least-squares deconvolution of apparent resistivity pseudosections. Geophysics 60, 1682 – 1690. Mufti, I.R., 1976. Finite-difference resistivity modeling for arbitrarily shaped two-dimensional structures. Geophysics 41, 62 – 78. Okabe, M., 1981. Boundary element method for the arbitrary inhomogenities problem in electrical prospecting. Geophys. Prospect. 29, 39 – 59.
292
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
Olayinka, A.I., Weller, A., 1997. The inversion of geoelectric data for hydrogeological applications in crystalline basement areas of Nigeria. J. Appl. Geophys. 37, 103 – 115. Peschel, G., 1967. A new favourable combination of resistivity sounding and profiling in archaeological surveying. Prospez. Archeol. 2, 103 – 115. Schleifer, N., Weller, A., Schneider, S., Junge, A., 2002. Investigation of a Bronze Age plankway by spectral induced polarization. Archaeol. Prospect. 9, 243 – 253. Schulz, R., 1985. The method of integral equation in the direct current resistivity method and its accuracy. J. Geophys. 56, 192 – 200. Schwarz, H.R., 1991. Methode der finiten Elemente: eine Einfqhrung unter besonderer Berqcksichtigung der Rechenpraxis, 3rd ed. Teubner, Stuttgart. Tsourlos, P.I., Szymanski, J.E., Tsokas, G.N., 1999. The effect of terrain topography on commonly used resistivity arrays. Geophysics 64, 1357 – 1363.
Tuyen, D., Canh, T., Weller, A., 2000. Geophysical investigations of river dikes in Vietnam. Eur. J. Environ. Eng. Geophys. 4, 195 – 206. Ullrich, B., Meyer, C., 2003. Geophysical prospection of dikes. Proceedings of Vietnamese–German workshop on Dike Monitoring, Hanoi, pp. 25 – 30. Ward, S.H., 1990. Resistivity and induced polarization methods. In: Ward, S.H. (Ed.), Geotechnical and Environmental Geophysics, Review and Tutorial, vol. I. Society of Exploration Geophysicists, Tulsa, OK, pp. 147 – 189. Weller, A., Seichter, M., Kampke, A., 1996. Induced-polarization modelling using complex electrical conductivities. Geophys. J. Int. 127, 387 – 398.
The effect of dike geometry on different resistivity configurations Thomas Henniga,T, Andreas Wellera, Tran Canhb a
Institut fu¨r Geophysik, Technische Universita¨t Clausthal, Arnold-Sommerfeld-Str. 1, D-38678 Clausthal-Zellerfeld, Germany b Institute of Geological Sciences, Vietnamese Academy of Science and Technology, Hanoi, Vietnam Received 19 March 2004; accepted 1 March 2005
Abstract Geoelectrical profiling with multi-electrode systems has become an important tool for monitoring dike embankments bordering rivers. Profiles running perpendicular to the dike axis are affected by the dike topography, with the amplitude of this effect dependent on the surface geometry and the choice of the electrode configuration. Investigations using seven different electrode configurations have shown that some configurations are less sensitive to the topography than others. The topography correction method (TCM) is an important tool for processing data from measurements at river dikes. This method is generally recommended for flank angles steeper than 108. The topography effect is calculated by two-dimensional finite element modelling. The resulting synthetic data of a homogeneous dike body are used to apply a topographic correction for each measurement. The topographic effect and correction procedure is demonstrated for synthetic dike data and for a data set from a river dike in Thai Binh province (Vietnam). The topography can be ignored for flank angles less than 258 if an averaged Half-Wenner electrode configuration is used. This configuration has proved to be less affected by undulated topography and the focusing effect of averaging the two data sets provides reliable structural information without the need for time-consuming data inversion. D 2005 Elsevier B.V. All rights reserved. Keywords: Electrical resistivity; Finite element method; Topography correction; Dike monitoring
1. Introduction Multi-electrode arrays are widely used in monitoring dike embankments that border rivers to prevent flooding (Ullrich and Meyer, 2003). Cables
T Corresponding author. Fax: +49 5323 722320. E-mail address: [email protected] (T. Hennig). 0926-9851/$ - see front matter D 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.jappgeo.2005.03.001
connect a large number of equally spaced electrodes with a control unit. The active electrodes are selected by the control unit according to a predefined measurement sequence. A set of different electrode configurations which are commonly used in multi-electrode surveys is shown in Fig. 1 with A and B being the current electrodes, M and N the potential electrodes, a the fixed electrode distance and n an integer factor. Each configuration has its
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
279
(a) Dipole - Dipole B
A
M
a
N
na
a
(b) Wenner-α A
M
N
na
na
B na
(c) Wenner-β B
A
M
na
na
N na
(d) averaged Wenner-α/β A
M
N na
na B
A
B na
M
N
(e) Half-Wenner forward A
M na
N
B
oo
A
B
oo
B
oo
na
(f) Half-Wenner backward N
M na
na
(g) averaged Half-Wenner A
M
N
na N
na M
A
(h) Schlumberger A
M na
N a
B na
Fig. 1. Commonly used electrode configurations for a multi-electrode survey. (a) Dipole–dipole. (b) Wenner-a. (c) Wenner-h. (d) Averaged Wenner a/h. (e) Half-Wenner forward. (f) Half-Wenner backward. (g) Averaged Half-Wenner. (h) Schlumberger.
advantages and drawbacks (Ward, 1990). The Wenner-a (Fig. 1b) configuration is suitable to resolve resistivity changes with depth. Gradient arrays like dipole–dipole configuration (Fig. 1a) provide a better lateral resolution but they become
severely affected by noise as the electrode spacing increases. Earlier investigations have shown that a special type of three electrode configuration which is called Half-Wenner configuration (Fig. 1e,f) provides some
280
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
advantages for resistivity profiling with multi-electrode systems (Peschel, 1967): (i) The number of measurements with a fixed number of electrodes is larger than for Wenner-a and dipole–dipole configuration. Since more levels can be measured an increased depth of penetration can be reached. (ii) For the infinite electrode B a good location for current injection can be found a large way outside the profile. (iii) Two Half-Wenner readings can be combined to a Wenner-a reading. (iv) The simple combination of a Half-Wenner forward configuration and a Half-Wenner backward configuration with a fixed position of electrode M has proved to act as a focussing tool (Kampke et al., 1998). The average of both readings, which is called the averaged HalfWenner (Fig. 1g), results in a pseudo-section that shows the main features of the subsurface resistivity distribution. It will be shown in this paper, that the averaged Half-Wenner configuration is also less sensitive than other commonly used configurations (e.g. Wenner-a, dipole–dipole) to a rough topography. Resistivity profiling is normally applied along the dike axis. Profiles are measured either at the top of the dike, the flanks, or even at the foot of the dike to explore the foundation. If some resistivity anomalies are found it is often difficult to predict the exact location of the anomalous structure. A two-dimensional interpretation assumes that the anomalies are caused by structures directly beneath the profile. But it is known that structures at both sides of the profile also affect the measurements. In order to find out the exact location of the zone of anomalous resistivity it is advisable to measure a perpendicular profile, which provides a cross-section of the dike body. Such a profile will usually be severely affected by the dike geometry. Different authors have investigated the influence of topography on resistivity readings (e. g. Fox et al., 1980; Tsourlos et al., 1999). Fox et al. (1980) have shown that a dipole–dipole survey crossing a ridge of homogeneous resistivity causes an apparent resistivity maximum under the top and minima at the two feet of the ridge. The maximum exceeds the intrinsic
resistivity by the factor 2.6. The minimum drops to 50% of the true resistivity. Tsourlos et al. (1999) used a topography similar to the dike geometry. The resulting resistivity pseudo-section for a dipole–dipole array shows a strong resistivity maximum under the dike body and adjacent minima near the feet of the dike. The values of the maximum and the minima are comparable to the model of Fox et al. (1980) that uses the same slopes of 308 for both flanks. Tsourlos et al. (1999) stated that slopes larger than 108 cause misleading artificial errors that have to be removed. In this paper, the effect of topography for different electrode configurations is demonstrated for both a model of synthetic data and field data from a Vietnamese dike. A topography correction method is applied (Fox et al., 1980) that largely eliminates the influence of the dike geometry from the measured data.
2. Modelling with finite element method (FEM) Forward modelling of geoelectrical measurements can be done by different numerical methods: boundary element method (e.g. Dieter et al., 1969; Okabe, 1981; Schulz, 1985), finite difference modelling (e.g. Mufti, 1976; Dey and Morrison, 1979a,b; Weller et al., 1996) or finite element method (FEM). Coggon (1971) was one of the first authors to use FEM for two-dimensional problems. An overview can be found in Hohmann (1988). The FEM provides some advantages for modelling geoelectrical data for an environment with a rough topography as the triangular elements are more flexible for more complicated structures compared with the rectangular elements of finite difference modelling. The Neumann-type boundary condition at the earth’s surface is automatically fulfilled by the formulation of the variation problem. Additionally, the discretisation errors at the surface are smaller than the errors using finite difference modelling due to the correct implementation of the boundary plain between airspace and subsurface for undulated topography. The actual problem of geoelectrical modelling is three-dimensional (3-D). The electrical potential U(x,y,z) in a non-uniform isotropic medium is governed by the following differential equation j½rð x; y; zÞU ð x; y; zÞ þ jjs ð x; y; zÞ ¼ 0;
ð1Þ
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
where the function r(x,y,z) describes the conductivity distribution and js(x,y,z) the source distribution of current density. A Neumann-type boundary condition at the surface and asymptotic conditions at the other model boundaries complete the problem. The differential Eq. (1) is valid for steady current flow and is also acceptable for alternating current at low frequencies if electromagnetic effects can be ignored. Dike structures can be represented by 2-D models. If the y-axis of the model is chosen to be parallel to the strike direction, the conductivity distribution becomes independent of y and can be written as r(x,y,z) = r(x,z). Even a 2-D resistivity structure causes a 3-D potential distribution. A pure 2-D modelling requires line sources in the strike direction (Mufti, 1976). However, most of the current electrodes used in geoelectrical prospecting can be assumed to be point sources, and the problem remains 3-D. The numerical expense in solving such problems can be reduced by using a Fourier-cosine transformation of potentials in strike direction y (Dey and Morrison, 1979a): Z l ˜ U ð x; k; zÞ ¼ 2 U ð x; y; zÞcosðkyÞdy; ð2Þ 0
where k denotes the wavenumber. The current source is assumed to be located in the plane y = 0. The application of the Fourier-cosine transformation to the 3-D Eq. (1) results in a 2-D partial differential equation: B B B B r U˜ þ r U˜ Bx Bx Bz Bz ð3Þ 2 ˜ rk U þ Idð x xs Þdð z zs Þ ¼ 0; where I is the current which is injected at position (x s,0,z s) and d is the Dirac delta function. Eq. (3) has to be solved for a set of discrete wavenumbers k. The associated boundary condition for the upper boundary is a Neumann-type condition with r air = 0. A mixed boundary condition as proposed by Dey and Morrison (1979b) is applied at the left, right and lower borderline: rndjU bU ¼ 0;
ð4Þ
with n being the unit vector perpendicular to the boundary and b a factor that defines the asymptotic behaviour of the potential field in a homogeneous
281
medium. Compared with the often used Dirichlet condition, it allows a higher accuracy with a smaller modelling area (Weller et al., 1996). Consequently, fewer elements and less calculating time are required. The discretisation of the subsurface is shown in Fig. 2. This model represents an idealized dike with flank angles of 308. The highest density of nodes is chosen in the inner area around the electrodes. The modelling area is extended to the left and right and to the bottom to satisfy the mixed boundary conditions. The inter-nodal distance increases continuously outside the inner area. Each quadrangular element is cut into four triangles. This implements an additional node in the centre of the element and consequently a higher level of discretisation of the model area. Using the static condensation procedure (Schwarz, 1991), the additional nodes are eliminated but the higher level of discretisation is kept. Thus, a better discretisation can be attained without increasing the dimension of the system of equations, while the calculation time is reduced. ˜ (x,k,z) is assumed to The transformed potential U vary linearly over each of the triangular elements. Combining the resulting equations for all elements yields a symmetric system of equations with a dimension equal to the number of nodes. A characteristic property of FEM is that the resulting matrix is symmetric, sparse, banded and diagonally dominant. The structure depends on the order in which the nodes are numbered. The bandwidth of the resulting matrix, which is equal to the maximal index difference of nodes in all elements, should be minimized. The properties of the matrix are used to solve the linear system of equations with the Cholesky algorithm (Schwarz, 1991) for all wavenumbers. The resulting ˜ (x,k,z) in the wavenumber domain are potentials U used to reconstruct the potential U(x,y,z) by an inverse Fourier-cosine transform. If the potential electrodes are located in the plane y = 0 the inverse Fourier-cosine transformation is performed by a simple integration Z 1 l ˜ U ð x; y; zÞ ¼ U ð x; k; zÞdk; ð5Þ p 0 which is done numerically by a combination of Gaussian quadrature and Laguerre integration as proposed by LaBreque et al. (1996). The choice of the wavenumbers affects the accuracy and the computing time of the algorithm. In our algorithm, nine wave
282
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292 4
2
electrode position finite element grid node
z in m
0
-2
-4
-6
-8 -4
-2
0
2
4
x in m
6
8
10
12
14
Fig. 2. Ridge model for synthetic modelling after Fox et al. (1980) with the discretisation used for FE modelling.
numbers are used in order to keep the error of the inverse transformation less than 0.5% (Kemna, 2000). For each electrode configuration the potential difference between the electrodes M and N and the geometry factor K are used to calculate the apparent resistivity for each measurement: qa ¼
UM UN K: I
ð6Þ
The forward problem is solved for each current source separately. The apparent resistivities for all relevant electrode configurations can be computed easily by a superposition of several pole–pole configurations. This method reduces the number of forward problems to the number of current electrodes. The number of measured configurations is generally much higher.
3. Topographic correction method (TCM) Topographic effects occur if the surface is not flat. The intensity of the effect depends first of all
on the roughness of the surface and additionally on the electrode configuration used. Assuming a flat surface with a homogeneous subsurface, a current injected at an electrode at the surface causes a potential field around the electrode where the isopotential planes are represented by concentric halfspheres in the subsurface. According to the Neumann-type boundary condition no current can enter into the air (r air = 0). Consequently, the potential lines are perpendicular to the surface. Considering an uneven topography, the potential planes have to act in the same way, they remain perpendicular to the undulating surface. Consequently, the subsurface potential distribution is disturbed. If this effect is ignored anomalies in the subsurface would not be correctly interpreted. In order to determine the topographic effect of a certain configuration the potential distribution for any current electrode position is calculated, using a constant resistivity q mod, for all elements of the model. The apparent resistivity q cal results from the potential difference between the positions of electrodes M and N using the flat surface geometric factor
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
K that ignores the small changes in electrode distances caused by the topography. In the case of a flat surface for all configurations the apparent resistivity is expected to be close to the constant resistivity q mod attributed to all elements. Slight differences, that are called discretisation errors, are caused due to the finite accuracy of the finite element algorithm. Assuming an undulating surface, the resulting apparent resistivity q cal may differ considerably from the intrinsic value q mod. In order to correct the resulting pseudo-sections for the topographic effect correction factors CF (Fox et al., 1980) are determined by dividing the constant resistivity value q mod, that is attributed to all elements, by the calculated apparent resistivity q cal: CF ¼
qmod : qcal
ð7Þ
The CF that are calculated for a given surface topography can then be applied to the measured apparent resistivity data. Only if surface topography alone causes a resistivity anomaly the TCM yields exactly the data that would have been measured on a flat surface. Considering both inhomogeneous resistivity structures and an undulating surface, the topographic correction by the CF assumes that topography response and the signal resulting from the resistivity target are linearly superimposed. Since the resistivity forward problem is nonlinear this assumption can only be an approximation. But comprehensive tests with simple resistivity structures and surface geometry have shown that the TCM yields satisfactory results (Fox et al., 1980; Tsourlos et al., 1999).
4. Synthetic modelling experiments 4.1. Comparison with the ridge model The developed algorithm is tested in comparison with the results of Fox et al. (1980). Therefore, an almost identical ridge model is used to determine the apparent resistivity for a homogeneous subsurface. The grid used is shown in Fig. 2. It consists of 405 nodes and 370 elements. The calculated results match the results of Fox et al. (1980) for the dipole–dipole electrode configuration. Some deviations are found in
283
the shallow subsurface. This effect probably occurs because of the different discretisation of the subsurface air interface in both algorithms. Fox et al. (1980) use finite element technique but they fill up the modelling area to a horizontal line above the surface with elements of high resistivity values. Though Fox et al. (1980) do not exactly describe the dimension of the elements used, the average difference between the published values and our results reaches eight per cent. The largest differences (up to 20%) appear for the configurations with minimal spacing. Further tests with the dike model from Tsourlos et al. (1999) show a slightly better agreement with our modelling results. The average difference reaches 7% (Hennig, 2003). 4.2. Dike model with termite nest River dikes in the northern provinces of Vietnam have a total length of about 5000 km. The social and economic development of this region largely depends on the integrity of the dike system to prevent flooding. The safety of river dikes is affected by a diversity of problems. Sandy formations in the dike foundations cause seepage effects. Various termite species dig cavities for their nests in the dike body and decrease the stability of the dike (Tuyen et al., 2000). The efficiency of geophysical methods for dike monitoring has to be improved to detect such defects in the dikes. To study the problem of termite nests, an asymmetric dike model with a resistive cavity has been investigated. The dike model which is shown in Fig. 3 consists of a steeper flank with a slope of 308 on the left side and the other flank descends with an angle of 158. The cavity with a cross-section area of 4 m2 is placed 1.5 m below the top of the dike. It represents a typical termite nest with a resistivity value of 1000 V m (black box in Fig. 3). A resistivity of 25 V m is attributed to the surrounding silty dike material. The grid for FEM consists of 3941 nodes and 3756 elements. Two apparent resistivity data sets of a Wenner-a configuration have been generated by FEM modelling. The first data set corresponds to the homogeneous dike cross-section without a termite nest. The resulting pseudo-section which is displayed in Fig. 4a is shifted vertically according to the dike surface geometry. The depth of penetration is assumed to be
284
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
15
z in m
10
termite nest 1000 Ωm
5 dikebody 25 Ωm
0 -5 15
20
25
30
35
40
45
50
55
60
65
70
x in m Fig. 3. Dike model for synthetic modelling of a river dike with an embedded resistive structure (1000 V m) in the dike body (25 V m).
(a) z in m
10 5 0
15
20
25
30
35
40
45
50
55
60
65 31
x in m
30 29
(b)
28 10
z in m
27 26
5
24 0
23 22 15
20
25
30
35
40
45
50
55
60
65
21
x in m
20 19
(c) ρ in
z in m
10
m
5 0
15
20
25
30
35
40
45
50
55
60
65
x in m Fig. 4. Resistivity pseudo-sections of synthetic Wenner-a data of dike model in Fig. 3. (a) Model without termite nest, pseudo-section without topographic correction. (b) Model with termite nest, pseudo-section without topographic correction. (c) Model with termite nest, pseudo-section after application of topographic correction.
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
na/2 for all pseudo-sections. This is a good approximation for Wenner configurations (Barker, 1989). Though no nest was modelled, the pseudo-section shows a resistive structure at shallow depth beneath the top of the dike. But this anomaly and the other resistive and conductive features result only from the dike topography. The second data set was generated including the nest in the dike body. The resulting pseudo-section (Fig. 4b) shows a more accentuated resistive anomaly at the location of the nest. But it becomes obvious from the comparison of both pseudo-sections that a reliable conclusion whether a nest is located in the dike can not be drawn. The application of TCM to the second data set according to the previously described procedure yields the corrected pseudo-section which is shown in Fig. 4c. Though the lateral resolution of Wenner-a configuration is limited, an extended weak anomaly with resistivity values up to 28 V m becomes visible beneath the top of the dike. A similar anomaly would have been generated from a cavity below a flat surface. Most other anomalous structures which were related to the dike topography disappear after application of TCM. It should be noted that the application of TCM to the synthetic pseudo-section in Fig. 4a results in a homogeneous section with constant resistivity of 25 V m. Another data set that consists of both forward and backward Half-Wenner data has been generated by FEM modelling. A 2-D inversion of this data set was performed with the program DC2DSIRT (Kampke, 1996) which is based on a finite difference forward procedure (Weller et al., 1996) and a simultaneous iterative reconstruction technique. This robust inversion technique, which uses a back-projection as starting model, yields smooth models accentuating the main features of subsurface resistivity distribution after few iterations. The software which enables the use of arbitrarily chosen configurations has been successfully applied for hydrogeological (Olayinka and Weller, 1997) and archaeological (Kampke, 1999; Schleifer et al., 2002) projects. Several tests have proved that the results are comparable with those of commercial programs (e.g. Loke and Barker, 1995). Fig. 5a represents the inversion result after ten iterations. A resistive anomaly can be found in the deeper subsurface with a maximum value of approximately 50 V m. Almost no indication can be detected
285
in the shallow subsurface region of the termite nest. The topographic effect hides the target and causes a misleading interpretation. The termite nest would not have been detected. Applying TCM to the modeled data set and inverting the corrected data set yields the cross-section shown in Fig. 5b. A resistive anomaly is situated exactly at the predetermined position. The deeper structures, which were caused by the dike topography, disappear completely. Resistivity values close to 25 V m are calculated for the dike material. The resistive anomaly that represents the termite nest reaches about 30 V m. Though the maximal value does not correspond to the true value of 1000 V m, the indicated structure shows the termite nest in the correct position. These examples show the influence of topography on geoelectrical data, demonstrating that the interpretation of uncorrected data sets can be wrong, especially for small anomalies like termite nests. The data quality in the theoretical model is ideal whereas field measurements are often affected by noise. Consequently, TCM becomes more important for field data. 4.3. Sensitivity to topography of different electrode configurations The topographic surveying results of a dike at the right side of Tra Li river in Thai Binh province have been selected to demonstrate the sensitivity of different electrode configurations to a typical dike geometry. The dike cross-section and the FEM grid are shown in Fig. 6. A homogeneous resistivity distribution is assumed in the whole dike body and the subsurface. The height of the dike is about 10 m and the length of the profile crossing the dike is about 35 m. The flank angles reach up to 238. The correction factors (CF) are calculated for seven different configurations: ! ! ! ! ! ! !
dipole–dipole (Fig. 1a), Wenner-a (Fig. 1b), Wenner-h (Fig. 1c), averaged Wenner-a/h (Fig. 1d), Half-Wenner forward (Fig. 1e), Half-Wenner backward (Fig. 1f), and averaged Half-Wenner (Fig. 1g).
286
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
(a) z in m
10 5 0
15
20
25
30
35
40
45
50
55
60
65
x in m ρ in
m 18
22
26
30
34
38
42
46
50
54
(b) z in m
10 5 0
15
20
25
30
35
40
45
50
55
60
65
x in m ρ in
m 24
25
26
27
28
29
30
31
Fig. 5. Inversion results of synthetic Half-Wenner forward and backward data of dike model in Fig. 3. (a) Without topographic correction. (b) With topographic correction.
A summary of all CF and the total number of data are compiled in Table 1. The resulting pseudo-sections of CF are embedded in the dike geometry plots of Figs. 7 and 8. The influence of topography on the
resistivity readings can be seen in all configurations. However, the intensity and spreading of CF-values varies over a wide range for the considered configurations. The commonly used dipole–dipole config-
electrode position
z in m
5
finite element grid node
0
-5
-10 -25
-20
-15
-10
-5
0
x in m
5
10
15
20
Fig. 6. Topographic profile of river dike in Thai Binh province with discretisation.
25
30
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
287
Table 1 Values of correction factors (CF) for different electrode configurations for the dike geometry shown in Fig. 6 Dipole–dipole Wenner-a Wenner-h Half-Wenner forward Half-Wenner backward Averaged Half-Wenner Wenner-a/h averaged
Number of data
Minima of CF
Maxima of CF
Mean of CF
S.D. of CF
904 664 664 800 800 1600 1328
0.69 0.84 0.65 0.68 0.72 0.90 0.93
1.43 1.30 1.27 1.24 1.41 1.10 1.10
1.01 1.01 0.98 0.98 1.01 1.00 1.00
0.13 0.08 0.14 0.13 0.16 0.05 0.04
uration (Fig. 7a) is highly affected by topography. The CF values range from 0.69 to 1.43. According to Eq. (7) a CF of 1.43 corresponds to 43% deviation of the
apparent resistivity q cal from the true resistivity q mod. Beside the significant topography effect dipole–dipole data are generally more disturbed by noise, especially
(a) z in m
5
0
-10
-5
0
5
10
15
-10
-5
0
5
10
15
-10
-5
0
5
10
15
-10
-5
0
5
10
15
(b) z in m
5
0
(c) z in m
5
0
(d) z in m
5
0
correction factor
0.6
0.8
x in m
0.95
1.1
1.3
1.5
1.7
Fig. 7. Correction factors for different electrode configurations. (a) Dipole–dipole. (b) Wenner-a. (c) Wenner-h. (d) Averaged Wenner-a/h.
288
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
(a) z in m
5
0
-10
-5
0
5
10
15
-10
-5
0
5
10
15
-10
-5
0
5
10
15
(b) z in m
5
0
(c) z in m
5
0
correction factor 0.6
0.8
x in m
0.95
1.1
1.3
1.5
1.7
Fig. 8. Correction factors for different electrode configurations. (a) Half-Wenner forward. (b) Half-Wenner backward. (c) Averaged HalfWenner.
for larger dipole distances n. Therefore, the depth of penetration is normally restricted. The spreading of CF-values is slightly lower for Wenner-h (Fig. 7c), and the Wenner-a configuration yields a further reduction, where the CF-values range from 0.84 to 1.30 (Fig. 7b). Similar CF variations were obtained from the HalfWenner forward configuration that is shown in Fig. 8a. The Half-Wenner backward CF vary between 0.72 and 1.41 (Fig. 8b). Both Half-Wenner forward and backward configurations are also sensitive to topography. On closer examination of both diagrams it becomes obvious that the minima of the forward mode coincide with the maxima of the backward mode and vice versa. Consequently, the average of both images (Fig. 8c) results in considerably smaller CF varying between 0.91 and 1.10. Thus, the sensitivity to topography decreases for the averaged configuration compared to the original forward and backward modes.
The averaged Wenner-a and Wenner-h (Fig. 7d) as well as the averaged Half-Wenner configuration (Fig. 8c) are less sensitive to rough topography. The CFvalues vary from 0.90 to 1.10 for averaged Wenner-a and Wenner-h. As shown in Table 1, not only the interval but also the standard deviation (S.D.) of the CF-values becomes significantly reduced for the averaged configurations. The standard deviation can be regarded as topography induced noise that is added to the measured signal. Since the suppression of noise is a key criterion of data acquisition, configurations with low standard deviation in CF-values should be preferred. As mentioned above, the number of measurements for the same number of electrodes is higher for HalfWenner configurations (1600) than for Wenner-a/h (1328) in this example. The information density of the subsurface becomes more than 15% higher. Regarding this result and the above-mentioned advantages, it is recommended to use the averaged
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
289
device (GeoSys Leipzig) in Half-Wenner forward and backward electrode configurations. The profile length is 33 m. A constant electrode spacing of 0.5 m was used along the profile. The maximum level of investigation is n = 16. The uncorrected Half-Wenner forward, Half-Wenner backward, Wenner-a and averaged Half-Wenner pseudo-sections are presented in Fig. 9. The Wenner-a data were generated by averaging two Half-Wenner readings with the same electrode spacing na and a common reference point in the centre between the electrodes M and N. It can be seen that the asymmetric configurations, like HalfWenner forward and backward, result in asymmetric
Half-Wenner configuration for geoelectrical data acquisition. The vertical and lateral resolution of this averaged configuration satisfies the requirements of dike investigations.
5. Field data The effect of topographic correction is demonstrated in Figs. 9 and 10 for a field data set from Tra Li river in Thai Binh province (Vietnam). The dike topography is the same as in the example shown in Fig. 6. The data was acquired with a GSM 150
(a) z in m
5
0 -10
-5
0
5
10
15
20
-10
-5
0
5
10
15
20
-10
-5
0
5
10
15
20
-10
-5
0
5
10
15
20
(b) z in m
5
0
(c) z in m
5
0
(d) z in m
5
0
x in m a
in
m 12
18
26
37
55
80
Fig. 9. Pseudo-sections of field data from Tra Li river dike. (a) Half-Wenner forward. (b) Half-Wenner backward. (c) Wenner-a. (d) Averaged Half-Wenner.
290
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
(a) z in m
5
0
-10
-5
0
-10
-5
0
5
10
15
20
5
10
15
20
z in m
(b)5
0
x in m in
m 12
18
26
37
55
80
Fig. 10. Inversion result of field data from Tra Li river dike. (a) Without TCM. (b) With TCM.
pseudo-sections. Some maxima of the forward measurement correspond with minima of the backward measurement and vice versa. As shown in the previous chapter, this topographic effect is widely reduced by an averaging overlay of Half-Wenner forward and backward pseudo-section. Compared with the Wenner-a pseudo-section, which is more distorted by topography, the averaged Half-Wenner image is smoother. The measured data set, consisting of 1148 HalfWenner forward and backward readings, was processed with the 2-D inversion program DC2DSIRT. In a first run the original data were inverted for a flat surface, whereas for the second run the topography correction was applied to the data. The resulting resistivity models are shown in Fig. 10a and b. Since only slight differences between the two images are visible it becomes obvious that a topography correction can be neglected if both Half-Wenner forward and backward data are available. The comparison with the averaged pseudo-section in Fig. 9d verifies that the inversion does not change the image significantly. The high imaging quality of the averaged Half-Wenner pseudo-section results from the focusing effect of this configuration (Kampke, 1999). All relevant features become visible in both the pseudo-section and the inverted model. A resistive anomaly is found in the shallow subsurface of the land side flank of the dike at
x = 2 m. As an excavation shortly after the survey has proved, this anomaly is caused by an extended termite nest. Since no excavations were allowed at the river side dike flank, the origin of the other resistive anomaly at x = 7 m could not be verified.
6. Conclusions Geoelectrical multi-electrode measurements are a suitable tool for investigating river dikes although topographic effects often complicate the interpretation of the measured data. To overcome the problems caused by a rough topography either topographic correction methods (TCM) or less sensitive electrode configurations should be used. TCM is based on a finite element forward modelling that considers the real dike geometry. The dike body is assumed to be homogeneous. The resulting section of correction factors (CF) represents the topographic effect at the measured data. The TCM uses CF to remove topographic effects from measured data. The interval of CF-values and their standard deviation enable a quantitative comparison of the sensitivity to topography of different electrode configurations. It has been shown that an appropriate combination of different electrode configurations reduces the
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
sensitivity of the generated pseudo-sections to the topography. This effect was demonstrated with the averaged Half-Wenner configuration and the average of the Wenner-a and Wenner-h configurations. The investigations have shown that the topographic effect becomes larger as the flank angle increases. The topographic correction can be neglected for both averaged configurations if the flank angles are less than 258. Since this condition is fulfilled at most dikes there is a clear advantage of using these configurations for dike monitoring. Other configurations that are more affected by topography like dipole–dipole and Wenner-a require a topographic correction even at lower flank angles. The use of configurations that are less affected by topography makes the interpretation also less sensitive to uncertainties of the geometric surveying data. The larger the topographic effects of the chosen configuration the more exact the surface topography has to be surveyed. Thus, the expense for geometric surveying can be reduced if appropriate configurations are applied. This conclusion remains valid even if modern inversion procedures which include an automatic topographic correction are used. Another advantage of the averaged configurations is the resemblance between the pseudo-section and the section obtained by inversion that results from the focusing effect of these combined configurations. Without any time-consuming inversion a reliable image of the subsurface structures can be generated in the field. Termite nests and other defects in the dike body can be easily detected in the pseudo-section plot. The operator can quickly provide the relevant information that would allow more detailed lines to be measured right away. Considering the low sensitivity to topography, the high image quality and the good data coverage, a survey with a combination of Half-Wenner forward and backward configurations is highly recommended for dike monitoring.
Acknowledgments The presented investigations were made in the joint scientific project bDevelopment of a geoelectrical technology to find termite nests in river dikesQ
291
between the Institute of Geological Sciences at the Vietnamese Academy of Sciences and Technology and the Institute of Geophysics at Technical University Clausthal (Germany) that was sponsored by the Volkswagen-Stiftung (2001–2004). We thank for the assistance provided by the authorities of Thai Binh province which enabled our field measurements at a dike section at Tra Li river.
References Barker, R.D., 1989. Depth of investigation of collinear symmetrical four-electrode arrays. Geophysics 54, 1031 – 1037. Coggon, J.H., 1971. Electromagnetic and electrical modeling by the finite element method. Geophysics 36, 132 – 155. Dey, A., Morrison, H.F., 1979a. Resistivity modelling for arbitrarily shaped two-dimensional structures. Geophys. Prospect. 27, 106 – 136. Dey, A., Morrison, H.F., 1979b. Resistivity modeling for arbitrarily shaped three-dimensional structures. Geophysics 36, 753 – 780. Dieter, K., Paterson, N., Grant, F., 1969. IP and resistivity type curves for three-dimensional bodies. Geophysics 34, 615 – 632. Fox, R.C., Hohmann, G.W., Killpack, T.J., Rijo, L., 1980. Topographic effects in resistivity and induced polarization surveys. Geophysics 45, 75 – 93. Hennig, T., 2003. Modellierung von Polarisationseffekten mit finiten Elementen. Diplomarbeit. Institut fqr Geophysik, Clausthal-Zellerfeld. Hohmann, G.W., 1988. Numerical modeling for electromagnetic methods in geophysics. In: Nabighian, M.N. (Ed.), Electromagnetic Methods in Applied Geophysics, vol. 1. Invest. Geophys. Soc. Exp. Geophys, Tulsa, OK, pp. 313 – 363. Kampke, A., 1996. Modellierung and Inversion von Daten der induzierten Polarisation fqr zweidimensionale Verteilungen der elektrischen Leitf7higkeit. Diplomarbeit. Institut fqr Geophysik und Meteorologie, Braunschweig. Kampke, A., 1999. Focused imaging of electrical resistivity data in archaeological prospecting. J. Appl. Geophys. 41, 215 – 227. Kampke, A., Weller, A., Peschel, G., 1998. Focussing effect on an averaged three-electrode configuration. Proceedings of the IVth Meeting of EEGS-European Section, Barcelona, pp. 865 – 868. Kemna, A., 2000. Tomographic inversion of complex resistivity— theory and application. Bochum, Ruhr-Universit7t, Dissertation. LaBreque, D.J., Miletto, M., Ramirez, W., Owen, E., 1996. The effects of noise on Occam’s inversion of resistivity tomography data. Geophysics 61, 538 – 548. Loke, M.H., Barker, R.D., 1995. Least-squares deconvolution of apparent resistivity pseudosections. Geophysics 60, 1682 – 1690. Mufti, I.R., 1976. Finite-difference resistivity modeling for arbitrarily shaped two-dimensional structures. Geophysics 41, 62 – 78. Okabe, M., 1981. Boundary element method for the arbitrary inhomogenities problem in electrical prospecting. Geophys. Prospect. 29, 39 – 59.
292
T. Hennig et al. / Journal of Applied Geophysics 57 (2005) 278–292
Olayinka, A.I., Weller, A., 1997. The inversion of geoelectric data for hydrogeological applications in crystalline basement areas of Nigeria. J. Appl. Geophys. 37, 103 – 115. Peschel, G., 1967. A new favourable combination of resistivity sounding and profiling in archaeological surveying. Prospez. Archeol. 2, 103 – 115. Schleifer, N., Weller, A., Schneider, S., Junge, A., 2002. Investigation of a Bronze Age plankway by spectral induced polarization. Archaeol. Prospect. 9, 243 – 253. Schulz, R., 1985. The method of integral equation in the direct current resistivity method and its accuracy. J. Geophys. 56, 192 – 200. Schwarz, H.R., 1991. Methode der finiten Elemente: eine Einfqhrung unter besonderer Berqcksichtigung der Rechenpraxis, 3rd ed. Teubner, Stuttgart. Tsourlos, P.I., Szymanski, J.E., Tsokas, G.N., 1999. The effect of terrain topography on commonly used resistivity arrays. Geophysics 64, 1357 – 1363.
Tuyen, D., Canh, T., Weller, A., 2000. Geophysical investigations of river dikes in Vietnam. Eur. J. Environ. Eng. Geophys. 4, 195 – 206. Ullrich, B., Meyer, C., 2003. Geophysical prospection of dikes. Proceedings of Vietnamese–German workshop on Dike Monitoring, Hanoi, pp. 25 – 30. Ward, S.H., 1990. Resistivity and induced polarization methods. In: Ward, S.H. (Ed.), Geotechnical and Environmental Geophysics, Review and Tutorial, vol. I. Society of Exploration Geophysicists, Tulsa, OK, pp. 147 – 189. Weller, A., Seichter, M., Kampke, A., 1996. Induced-polarization modelling using complex electrical conductivities. Geophys. J. Int. 127, 387 – 398.