- Email: [email protected]
Water content and porosity estimated from ground-penetrating radar and resistivity
Journal of Applied Geophysics 58 (2006) 99 – 111 www.elsevier.com/locate/jappgeo
Water content and porosity estimated from ground-penetrating radar and resistivity Anita Turesson Department of Geology, Earth Sciences Centre, Go¨teborg University, Box 460, SE-405 30 Go¨teborg, Sweden Received 1 September 2004; accepted 29 April 2005
Abstract Both ground-penetrating radar and the resistivity method have proven to be useful tools for exploring water content variations, since related parameters such as dielectric constant and the resistivity of rocks and sediments are highly dependent on moisture content. These methods were used independently to estimate volumetric water content in the unsaturated zone and porosity in the saturated zone in a 100-m sandy section. Two sample sites along the profile were also chosen for a shallow geophysical investigation and soil sampling, to enable the calibration and verification of the indirect geophysical methods. The grain distribution at these sites is dominated by medium-sized sand (0.25–0.5 mm). The water content was 6.9 vol.% and calculated porosities are 37% and 40% respectively. At each of these sites the mean water content values calculated from resistivity are within one percentage unit of measured water content while those calculated from ground-penetrating radar give higher values by as much as 2.9 percentage units. The water contents in the unsaturated zone in the section, estimated from resistivity and ground-penetrating radar, show very similar trends, although that deduced from ground-penetrating radar is generally somewhat larger, consistent with the results from the sample sites. The mean porosity values obtained from the two methods in the saturated zone are in good agreement. D 2005 Elsevier B.V. All rights reserved. Keywords: Ground-penetrating radar; Resistivity; Volumetric water content; Porosity
1. Introduction Soil water content and porosity are important variables in hydrological processes and are of primary interest in hydrogeological investigations. Ground-penetrating radar (GPR) has proven to be a
E-mail address: [email protected]. 0926-9851/$ - see front matter D 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.jappgeo.2005.04.004
promising technique for estimating water content in soil (Greaves et al., 1996; Van Overmeeren et al., 1997; Huisman et al., 2001). Using GPR in combination with mixing formulae, the water content can be estimated from dielectric constants, which are calculated from interval velocities of radar waves. The GPR method has successfully been applied to shallow (less than 50 m) geological surveys (Davis and Annan, 1989). The advantage of the method is
100
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
its vertical and lateral resolution in high resistivity environments. The main restriction of the method is the limited penetration in conductive materials, such as clays or soils with saline or contaminated pore water. Empirical relationships relating electrical resistivity to porosity have long been known (Archie, 1942). The water content in a clean sand formation can be estimated with Archie’s saturation equation. The applicability of geoelectric methods is not restricted to certain sedimentary sequences. However, the principle of equivalence, i.e. the thickness and resistivity of a unit can vary within certain limits and still give equivalent models (Parasnis, 1997), makes it difficult to estimate true subsurface resistivities. GPR and geoelectric methods have been used in combination to estimate water content and porosity in different ways. Dannowski and Yaramanci (1999) used GPR to constrain the geometry used in geoelectric inversion and compared the results from the two methods. Garambois et al. (2002) combined GPR and geoelectrics to estimate water content and water conductivity variations in the unsaturated zone. In contrast to previous investigations, this study aims to independently evaluate two methods, GPR and resistivity, used to assess water content (porosity in the saturated zone) according to the relationship of Topp et al. (1980) and Archie’s saturation equation (Ward, 1990) respectively. The Topp et al. (1980) equation, relating the dielectric constant to water content, was chosen among other mixing formulae (Mavko et al., 1999) because of its simplicity and good results (Greaves et al., 1996). Using standard techniques, in combination with well known empirical relationships, additional parameters can easily be obtained. This study uses geophysical data collected in a 100 m profile during two consecutive days in stable weather conditions. Two spots along the profile were chosen for small-scale shallow investigations to allow calibration and verification of the indirect geophysical methods. These investigations included GPR measurements, vertical electrical soundings, and soil sampling at a depth of 1 m. The measurements and sampling were made during two consecutive days at each site within a month of the profile survey. The data from each site were compared separately. The soil samples were analyzed to determine grain size distribution and gravimetric
water content, which was converted to volumetric water content. In addition, the P-wave seismic refraction was used to discriminate between the unsaturated and saturated zone, to support the results obtained from the GPR and resistivity. The objective of this paper is to compare and evaluate the variations in water content (porosity in the saturated zone) in sand estimated independently by two standard geophysical techniques, ground-penetrating radar and resistivity, in combination with the relationship of Topp et al. (1980) and Archie’s saturation equation (Ward, 1990).
2. Geological setting The test site is located at Veddige 70 km south of Go¨teborg on the Swedish west coast at an elevation of 15 m (Fig. 1). The post-glacial marine limit in the region is 65 m above sea level (Pa˚sse, 1986). During the overall post-glacial regression a small transgression (the Tapes transgression) occurred in certain parts of southwestern Sweden, reaching 17 m above sea level in the area (Pa˚sse, 1986). Large quantities of mostly glaciofluvial sediments, originating from a terminal moraine in the vicinity, were redeposited in an old channel. The area is underlain by wave-sorted sand and gravel which overlies the
Fig. 1. Map of southwestern Sweden showing the location of the site discussed in the text.
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111 Table 1 Grain size distributions presented as the mean of two samples at each site Sample site
Grain size distribution (mm)
Weightpercent
Measured h
Calculated porosity
30 m
Gravel (N2) Sand (0.071–2) Silt (b0.071) Gravel (N2) Sand (0.071–2) Silt (b0.071)
16 80 4 21 75 4
6.9%
39.9%
6.9%
37.3%
60 m
101
m. The information from these sample pits is summarized in Table 1.
3. Seismic refraction 3.1. Method
Measured volumetric water content (h) is the mean of three samples. A grain density of 2.65 g/cm3 was used in the porosity calculation. The samples were taken at about 1 m depth. See Figs. 3 and 7 for locations of sample sites.
glacial marine clay that in turn overlies sand and gravel (Pa˚sse, 1986). The surveyed profile is located on the redeposited sand with a maximum depth of 14 m. Soil samples were taken at a depth of 1 m at 30 and 60 m along the profile. These soil samples were analyzed for gravimetric water content and converted to volumetric water content using soil and pore water density estimates between 1.58– 1.68 and 1.00 g/cm3 respectively. The two sample sites gave the same result of 6.9 vol.%. The calculated porosities, using a grain density of 2.65 g/ cm3, are 37% and 40%. The pits revealed a sharp boundary between soil and sand at 0.33 m. The sand was analyzed to determine particle size distribution and is mostly of medium size (0.25–0.5 mm), 58% at sample site 30 m and 42% at 60
To determine the depth to the water table by a third independent method, we made a P-wave refraction survey using an ABEM Terraloc Mark 6 seismograph. The data were collected using 36 geophones (10 Hz) with 2-m spacing. Two spreads were measured, making a total length of 142 m. The energy source was a sledgehammer hitting a steel plate and three blows with the hammer were stacked in each record. The data were processed using the delay-time method (Pakiser and Black, 1957) followed by ray-tracing (Yacoub et al., 1970). 3.2. Result The seismic model shows the interface between unsaturated and saturated zones declining from 8.1 to 10.1 m in the investigated area (outlined in Fig. 2). The first layer with 375 m/s velocity is dry sand. The second layer with 1465 m/s velocity is interpreted as wet sand. The P-wave velocity increases greatly when water saturation reaches 100% (Bachrach and Nur, 1998), the refraction method is therefore suitable to determine the water table depth in coarser material such as sand. The water table deduced from seismic
Fig. 2. The seismic refraction model showing the interface between unsaturated and saturated zones. The water table declines from 8.1 to 10.1 m in the outlined area which coincides with the GPR profile.
102
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
refraction is plotted in all further models calculated from GPR and CVES.
4. Ground-penetrating radar 4.1. Method Two different GPR techniques were used in this study, the common-offset and the common midpoint (CMP) methods. The former was used to get an overview of the subsurface and to reveal any steeply declining surface, while the latter was used to get velocities with depths at regular intervals along the profile. Additional CMP soundings were measured at the two sample sites at 30 and 60 m to enable comparisons between water contents deduced from GPR and known water contents. GPR data is usually collected with the commonoffset method, in which the transmitting and receiving antennae are a fixed distance apart (Davis and Annan, 1989). The 100 m common-offset profile was measured using a Ramac GPR system from Mala˚ Geoscience. Antennae were 2 m apart with a centrefrequency of 50 MHz. The antennae were orientated parallel to each other and perpendicular to the profile. The data were collected every 0.2 m at a sample frequency of 500 MHz. To improve the signal-tonoise ratio, every trace was vertically stacked using the auto-stack option, which means that every trace was stacked between 32 and 64 times. A bhip chainQ, calibrated to actual length, was used to measure distance. This is a simple way to control the distance between each trace measured. A cotton-thread runs out as the operator slowly walks the profile. Every 0.2 m (in this case) the radar is triggered and a radar pulse is transmitted and subsequently received by the antenna. A more labour-intensive way to obtain data is to make an entire GPR survey using the multi-offset CMP technique (Fisher et al., 1992). In this case the radar data can be used for traditional seismic processing and are sorted in CMP-gathers, which are used in velocity analyses. The one-dimensional velocity models analyzed from CMP-gathers are then interpolated for the construction of a two-dimensional velocity profile. However, a less laborious way, used in this survey, to obtain such a 2-D velocity profile is to conduct individual multi-offset CMP soundings at
appropriate intervals along the profile. These measurements can easily be collected in the field using the hip chain measuring device. In a multi-offset CMP sounding the separation between transmitting and receiving antenna is continually increased from a fixed central location while the two-way travel time to reflectors are measured. Any subsurface contrast in electromagnetic properties results in energy being reflected back to the surface. Each reflection measured in this manner is used to derive the RMS (root mean square) velocity down to it. The CMP soundings were conducted using 50 MHz antennae with a sampling frequency of 500 MHz. At each CMP location the antenna separation was increased from 0 to 20 m, with increments of 0.2 m. The measuring device was placed at the midpoint. The true distance walked was corrected in the processing of data. To ensure that the antennae were moved equally from the midpoint, the measurements were performed stepwise using a measuring tape to control distances. Two persons are needed for this procedure. The spacing of individual CMP soundings, which should be measured in a profile, is a compromise between lateral and vertical variations of radar wave velocity and the time and effort to make the measurements. In this case we chose to conduct a spacing of 5 m and a total number of 21 CMP’s, which took about 2 h to collect. At the sample sites, CMP soundings were conducted in two directions to reveal possible three-dimensional geometry. The antennae centre frequency used was 200 MHz with a sampling frequency of 2000 MHz. The trace increment was 0.1 m. The vertical auto-stack function was used in all CMP measurements. The processing of GPR data (50 MHz) included time-zero adjustments and low-cut filtering (dewow), which removes low-frequency induction effects on the radar equipment. The data were also compensated for geometrical spreading and attenuation. The linear part of the gain was set to 0.02 (1/pulse width) and the exponential part was set to 0.03 dB/m (Davies and Annan, 1989). AGC scaling was used for display. The semblance approach (Yilmas, 1987) was used to pick preliminary RMS (or normal move-out) velocities. If the CMP data contains many and closely-spaced (in time) reflections it could be difficult to distinguish between real reflections arising from the interface between two electrically different media, and just a
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
complicated reflected wavelet. The compensation for geometrical spreading and attenuation helped to pick strong reflections in true amplitude display. To refine the velocity picks, hyperbolae were superimposed on the actual CMP gather to attain optimal fit. 4.2. Water content deduced from GPR From CMP-gathers the RMS velocity to reflectors is determined. The interval velocity, between reflectors, is calculated using the Dix (1955) equation: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v2RMS;n tn v2RMS;n1 tn1 Vi ¼ ð1Þ tn tn1 Where Vi is the interval velocity, v RMS are the RMS velocity and t n and t n 1 is the reflected ray two-way travel times to the nth and (n1)th reflectors respectively. The GPR produces high frequency electromagnetic energy in the 10 to 1000 MHz range. The term used to describe the high frequency electromagnetic properties of materials is called the dielectric constant or the relative permittivity. The complex dielectric constant is given by: KT ¼ KV þ jfKW þ ðrdc =xe0 Þg
ð2Þ
Where KV is the real part of the dielectric constant, and KW is the imaginary part of the dielectric constant or the electric loss, r dc is the zero-frequency conductivity, x is the angular frequency, e 0 is the free-space permittivity (8.854 10 12 F/m), and j is ( 1)1 / 2. In soils where the electric loss is small, that is with a lower conductivity than 10 mS/m, the following
103
relationship can be used to calculate the real part of the dielectric constant (KV; Davis and Annan, 1989): c ð3Þ V c pffiffiffiffiffi KV Where V is the propagation velocity of an electromagnetic wave in a medium with a real dielectric constant of KV, c is the velocity of an electromagnetic wave in free space (3 108 m/s). Topp et al. (1980) found that the real part of the dielectric constant (KV) seems to be highly sensitive to volumetric water content, but only weakly sensitive to soil type and density. They used a wide range of soil samples, sandy loam to clay, to derive an empirical relationship between the apparent (measured) dielectric constant and volumetric water content: hv ¼ 5:3 102 þ 2:92 102 Ka 5:5 104 Ka2 þ 4:3 106 Ka3
ð4Þ
Where h v is the volumetric water content (the ratio of water volume to total sample volume). For low-loss materials K a c KV where K a is the apparent dielectric constant. The water content (h) equals the product of porosity (/) and water saturation (S w). In water saturated soils the water content (h) is a measure of porosity (/). h ¼ /dSw
ð5Þ
4.3. Results The 100 m long, GPR common-offset profile is presented in a time–depth section in Fig. 3. The
Fig. 3. The GPR common-offset profile with depth represented by time. The arrows show the water table between 20 and 100 m at about 150 ns and the lowest picked reflector used in velocity analyses between 215 and 260 ns. Sample locations are indicated.
104
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
Fig. 4. Examples of CMP measurements used in the velocity analyses with a one-dimensional velocity model (continuous line: interval velocity, dotted line: RMS velocity), hyperbolic adaptions and semblance images at a) 30 m, b) 60 m and c) 90 m along the profile.
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
105
Fig. 4 (continued).
distinct reflector between 20–100 m at about 150 ns is interpreted as the water table in sand, which also is supported by the seismic refraction model (Fig. 2). The slightly declining reflector below that, from 215 to 260 ns, marks the lowest reflector picked in velocity analyses. Note the diffractions from cables (known
locations) in the upper part at about 72 and 85 m. Examples of CMP gathers with semblance and calculated interval velocities are shown in Fig. 4. The 21 CMP gathers were used to construct the two-dimensional velocity section converted to depths (Fig. 5). This section shows interval velocities between 72 and
Fig. 5. The GPR two-dimensional interval-velocity profile constructed from 21 one-dimensional velocity models. The continuous line shows the water table deduced from seismic refraction.
106
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
Fig. 6. Volumetric water content section calculated from GPR interval velocities using the relationship of Topp et al. (1980). The continuous line shows the water level deduced from seismic refraction.
149 m/As. The line marks the water table deduced from the seismic method. The mean interval velocity in the saturated zone is 76 m/As. The dielectric constants were calculated using Eq. (3) and the Topp et al. (1980) relationship (Eq. (4)) was used to convert the dielectric constants to volumetric water content (Fig. 6). The main part of the unsaturated zone shows water content between 10% and 12% but there are areas with lower water content (7–10%). The western part deviates with higher water content. In the saturated zone the water content varies between 22% and 31% (mean 28%). The CMP soundings at the sample sites revealed no three-dimensional heterogeneity. The interval velocities 127 m/As (at 30 m) and 125 m/As (at 60 m) yielded water contents of 9.4 and 9.8 vol.% respectively. These water contents differ by up to 2.9 percentage units from measured water content at these sites (Table 2).
5. Resistivity 5.1. Method The resistivity method is based on measuring the electrical potential which results from an applied direct electrical current flowing in the ground. The distribution of the electrical potential field depends in turn on the resistivity of the ground. In a multi-electrode array the measured result is displayed as a two-dimensional variation of apparent resistivity. Software applications for inversion of two-dimensional apparent resistivity to solve for true resistivity can be classified as either smooth inversion (DeGroot-Hedlin and Constable, 1990) or block inversion (Inman, 1975) methods, each of which has some disadvantages. Smooth inversion has a tendency to smear both resistivity and depth to interfaces even in the case of well-defined structures with sharp resistivity contrasts. On the other hand
Table 2 GPR wave velocity and vertical electrical soundings (VES) measured north–south (N–S) and east–west (E–W) directions Sample site
GPR velocity m/As
Calculated h (%)
VES q (Vm)
RMS error (%)
Calculated h (%)
30 30 60 60
127 127 125 125
9.4 9.4 9.8 9.8
11,300–11,500 –11,700 7500–7600–7800 6500–6600–6800 5600–5800–6100
0.9 1.3 1.6 3.2
5.3–5.3–5.2 6.8–6.7–6.6 7.4–7.3–7.2 8.1–7.9–7.7
m (N–S) m (E–W) m (N–S) m (E–W)
Equivalent VES models resistivities (q) up to 1.2% fit, and the RMS (root mean square) errors for the best model (middle value) are given. To calculate the volumetric water content (h) from Archie’s formula (Eq. (7)), pore water conductivity of 13 mS/m and m = n = 1.7 was used. See Figs. 3 and 7 for locations of sample sites.
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
block inversion requires a starting model close to the truth, which is rarely known precisely and is difficult to construct especially in complex cases. We decided to use the smooth inversion routine. The continuous vertical electrical soundings (CVES) were conducted with an ABEM Lund imaging system, using Wenner measurement geometry. The electrode spacing varied from 3 to 72 m in the 240 m profile. The CVES profile was centred on the GPR profile. A software package based on the smoothnessconstrained least-square method (DeGroot-Hedlin and Constable, 1990) was used to invert apparent resistivity to true resistivity. At the sample sites, 30 and 60 m, vertical electrical soundings (Schlumberger geometry) were measured in two directions, using the ABEM equipment. A total of 13 measurements were done moving the current electrodes equally and stepwise from the midpoint, with distances increasing from 0.5 to 8 m between the midpoint and current electrodes. 5.2. Water content deduced from resistivity Electrical conduction in soil is largely electrolytic, taking place in connected pore spaces and along grain boundaries. The relationship between resistivity (the inverse of conductivity) and porosity in sedimentary clay-free rocks is expressed by the formation factor, which is the ratio of the resistivity of the porous media to that of the pore fluid (Archie, 1942; Ward, 1990). q F¼ ¼ ad/m ð6Þ qw Where F is the formation factor, q is the bulk resistivity of the rock, q w is the resistivity of the
107
pore fluid, / is the porosity, and a and m are constants. A general form of Archie’s saturation equation is: q ¼ qw /m Swn
ð7Þ
Where q and q w is the bulk resistivity of the rock and the resistivity of the water respectively, / is the porosity, S w is the fractional water saturation, and n is the saturation exponent, which normally is equal to 2 (Ward, 1990). For a water-saturated rock Eq. (7) is reduced to: q ¼ qw /m
ð8Þ
Jackson et al. (1978) found that the exponent m was dependent on the shape of the particles, increasing as they became less spherical, while variation in size appeared to have little effect. Samples of natural sand have values of m in the range 1.4 to 1.6. In solving Archie’s formula (Eq. (7)) for volumetric water content, which is the product of porosity and water saturation, parameters such as the electrical conductivity of pore water, m and n had to be estimated. To do this, volumetric water content was calculated using a range of conductivities (2–40 mS/ m) and m = n(1.3–2). After comparison between calculated and measured water content from the two sample sites conductivity was estimated at 13 mS/m and m = n = 1.7. This conductivity for pore water is supported by measurement in a nearby well which was also 13 mS/m. In this investigation it was assumed that these parameters would not change in the section.
Fig. 7. The resistivity model. The outlined area coincides with the GPR profile and is used in the water content calculation. Sample locations are indicated.
108
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
5.3. Results
6. Discussion and conclusion
The 240 m long resistivity profile is presented in Fig. 7. The area investigated and compared to GPR measurements is outlined. Two anomalies with lower resistivities are seen at about 72 and 85 m. These are probably due to artefacts from cables, also seen in the GPR section (Fig. 3). Although the cables are less than 1 m below surface, the effects of resistivity distortion are seen to considerably greater depths. The water content was calculated using Archie’s equation (Eq. (7)), with water resistivity set to 77 Vm (the inverse of 13 mS/m) and m = n = 1.7. The resistivitybased water content section (h resistivity; Fig. 8) is more varied than the corresponding GPR section (h GPR; Fig. 6), containing both lower and higher values, although the general trend is the same with the exception of the two artefacts mentioned above. The water content in the unsaturated zone is between 5% and 14% with the lowest value about 2 percentage units lower than in corresponding h GPR profile. The increasing water content in the western part can also be seen in the h GPR profile. In the saturated zone the water content varies between 16% and 40% (neglecting the anomalies). The electrical soundings at the sample sites show some anisotropy, especially at 30 m. At each site mean water content values calculated from resistivity, 6.0% (30 m) and 7.6% (60 m), are within 1 percentage unit of the measured water content, 6.9% (Table 2).
Referring to data from the sample pits (Tables 1 and 2) the volumetric water content calculated from the Topp et al. (1980) equation is higher than both the measured value and that calculated from the Archie equation. Considering the corresponding profiles, h GPR and h resistivity (Figs. 6 and 8), the same trend is true for the larger part of the section, which can be seen in Fig. 9 showing the difference between these two results. However, except for the two anomalies the h resistivity shows higher water content near to the surface. This could be due to the fact that the smallest electrode spacing is 3 m which enables a shallower measurement of water content than for h GPR where the first picked reflector in the velocity analysis is at a depth of about 4 m. The h resistivity profile seems to be more detailed due to a denser sample grid with depth, but also probably due to the smooth inversion routine, which renders a gradual change of resistivities even if sharp boundaries exist in the subsurface (Olayinka and Yaramanci, 2000, 2002). On the other hand the less detailed h GPR is restricted by the number of reflections present in the subsurface and by the distance between the CMP’s. So in general, the difference in data density accounts for some of the differences between the two methods seen in Fig. 9. In the saturated zone the water content, which is the porosity when the pores are saturated, varies
Fig. 8. The volumetric water content section calculated from resistivity using Archie’s (1942) equation. Water conductivity is set to 13 mS/m and m = n = 1.7. The continuous line shows the water table deduced from seismic refraction.
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
109
Fig. 9. This section shows the difference between volumetric water contents deduced from GPR and resistivity.
between 22% and 31% in h GPR (Fig. 6) with a mean of 28%. The corresponding range for h resistivity (Fig. 8) is wider, roughly between 16% and 40% but also with a mean of 28%, calculated between depths of 10.80 and 12.89 m (excluding the anomaly between 82 and 90 m). In the sample pits the porosity was measured at 37% and 40% respectively. This suggests a decrease in porosity by approximately 10 percentage units at about 10 m depth. However, the measured porosity from 1 m depth is the total porosity which is the sum of effective, trapped and isolated porosities whereas the porosity calculated at 10 m depth from the resistivity method is the effective porosity because the current is largely electrolytic. The GPR method is based on wave propagation by analogy with seismic methods, so in that sense the porosity calculated from GPR is the total porosity. As the methods give very similar results this would indicate that the effective porosity equals the total porosity in this case and that no trapped and isolated pores exist. Porosity is governed by many factors such as the uniformity of grain size, (sorting), grain shape, packing, and compaction during and after deposition. Packing alone can contribute significantly to the difference in porosity. The end members of packing modes for spheres of uniform size, the cubic and rhombohedral packing, have porosities of 48% and 26% respectively (Graton and Fraser, 1935). Sorting also has a large influence on the porosity: up to 25% difference between well-sorted and very poorlysorted sands of the same mean grain size were reported by Beard and Weyl (1973). The dominant
factors in this study are not known but if the results indicated by two methods are correct, it is likely that more than one factor are responsible for the relatively large decrease in porosity with depth. Unfortunately, no direct control is available to verify the porosity at 10-m depth. A precondition for using the GPR in this type of study is the presence of several electromagnetically contrasting layers. This is not often a problem as sandy material is commonly stratified. Another requirement is that the layers are more or less horizontal since the equation of Dix (1955) is valid only for horizontal surfaces. The Topp et al. (1980) equation requires values for dielectric constant for each layer, derived from the interval velocity, as input for water content calculation. Thus the RMS velocity has to be carefully picked for each layer. In this study these were determined using both the semblance approach to select strong reflections and hyperbolae fitted to CMP gathers to refine the two-way times and RMS velocities. These velocities were thoroughly analysed and the maximum estimated error could be F1 ns in two-way-time and F0.003 m/ns in RMS velocity, which would result in about 1.5% error in water content. The smooth inversion method of DeGroot-Hedlin and Constable (1990), used for the geoelectrical data, has the advantage of being fully automatic and not needing any prior information, however the gradual change from high to low resistivities which is inherent to this method makes it difficult to determine an intrinsic value for porosity. A prerequisite for a reliable water-content model deduced
110
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
from resistivity is that true resistivities can be derived. This could also be problematic due to the equivalence principle relating depth and resistivity. An alternative inversion routine, the block inversion (Inman, 1975), has proven to be useful when the subsurface layer geometry is simple (Dannowski and Yaramanci, 1999). In more complex environments the input model required for this inversion routine may be too difficult to estimate to get a satisfying result. In using the Archie equation several parameters, such as the constants m and n and sometimes also water conductivity, have to be estimated. The small-scale investigation at the sample sites was used to calibrate these parameters. Furthermore, it was assumed that these parameters did not change in the section. In this survey this assumption can be justified by the facts that the constants m and n are related to pore shape and pore fluid and that the wave-sorted sandy material in the section has the same origin and genesis as sand washed out from a terminal moraine. In summary, two methods, ground-penetrating radar and resistivity, were independently evaluated for their capability to assess water content and porosity in a sandy section. The methods were used in combination with empirical relationships. Additional information from two sample sites was used to relate the indirect methods to known water content and porosity. This is of special importance when using empirical relationship in different specific environments. The results obtained showed very similar trends of water-content distribution, although absolute values differ somewhat, and there is a good agreement between the methods in the saturated zone if the mean porosity is compared. The use of two independent methods greatly strengthens the results which can be obtained in this type of study.
Acknowledgement I would like to thank associate professor Gustaf Lind, Earth Sciences Centre, Go¨teborg University, for assistance during fieldwork and for constructive review of the manuscript. I also thank professor David Cornell, Earth Sciences Centre, Go¨teborg University, for correction of the English language.
References Archie, G.E., 1942. The electrical resistivity log as an aid in determining some reservoir characteristics. Transactions of the American Institute of Mining and Metallurgical Engineers/Petroleum Division 146, 54 – 62. Bachrach, R., Nur, A., 1998. High-resolution shallow-seismic experiments in sand: Part I. Water table, fluid flow, and saturation. Geophysics 63 (4), 1225 – 1233. Beard, D.C., Weyl, P.K., 1973. The influence of texture on porosity and permeability on unconsolidated sand. The American Association of Petroleum Geologists Bulletin 57, 349 – 369. Dannowski, G., Yaramanci, U., 1999. Estimation of water content and porosity using combined radar and geoelectrical measurements. European Journal of Environmental and Engineeering Geophysics 4, 71 – 85. Davis, J.L., Annan, S.P., 1989. Ground-penetrating radar for highresolution mapping of soil and rock stratigraphy. Geophysical Prospecting 37 (5), 531 – 551. DeGroot-Hedlin, C., Constable, S., 1990. Occam’s inversion to generate smooth, two-dimensional models from magnetotelluric data. Geophysics 55 (12), 1613 – 1624. Dix, C.H., 1955. Seismic velocities from surface measurements. Geophysics 20 (1), 68 – 86. Fisher, E., McMechan, G.A., Annan, A.P., 1992. Acquisition and processing of wide-aperture ground-penetrating radar data. Geophysics 57 (3), 495 – 504. Garambois, S., Se´ne´chal, P., Perroud, H., 2002. On the use of combined geophysical methods to assess water content and water conductivity of near-surface formations. Journal of Hydrology 259, 32 – 48. Graton, L.C., Fraser, H.J., 1935. Systematic packing of spheres, with particular reference to porosity and permeability. Journal of Geology 43, 785 – 909. Greaves, J.G., Lesmes, D.P., Lee, M.L., Tokso¨z, M.N., 1996. Velocity variations and water content estimated from multi-offset, ground-penetrating radar. Geophysics 61 (3), 683 – 695. Huisman, J.A., Sperl, C., Bouten, W., Verstraten, J.M., 2001. Soil water content measurements at different scales: accuracy of time domain reflectometry and ground-penetrating radar. Journal of Hydrology 245, 48 – 58. Inman, J.R., 1975. Resistivity inversion with ridge regression. Geophysics 40 (5), 798 – 817. Jackson, P.D., Taylor, S.D., Stanford, P.N., 1978. Resistivity–porosity particle shape relationships for marine sands. Geophysics 43 (6), 1250 – 1268. Mavko, G., Mukerji, T., Dvorkin, J., 1999. The rock physics handbook. Tools for Seismic Analysis in Porous Media. Cambridge University Press. 329 pp. Olayinka, A.I., Yaramanci, U., 2000. Use of block inversion in the 2-D interpretation of apparent resistivity data and its comparison with the smooth inversion. Journal of Applied Geophysics 45, 63 – 81. Olayinka, A.I., Yaramanci, U., 2002. Smooth and sharp-boundary inversion of two-dimensional pseudosection data in presence of a decrease in resistivity with depth. European Journal of Environmental and Engineering Geophysics 7, 139 – 165.
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111 Pakiser, L.C., Black, R.A., 1957. Exploring for ancient channels with the refraction seismograph. Geophysics 22 (1), 32 – 47. Parasnis, D.S., 1997. Principles of Applied Geophysics. Chapman & Hall, UK. 429 pp. Pa˚sse, T., 1986. Beskrivning av jordartskartan Kungsbacka SO. Sveriges Geologiska Underso¨kning Ae 56 (106 pp.). Topp, G.C., Davis, J.L., Annan, A.P., 1980. Electromagnetic determination of soil water content: measurements in coaxial transmission lines. Water Resources Research 16 (3), 574 – 582. Van Overmeeren, R.A., Sariowan, S.V., Gehrels, J.C., 1997. Ground penetrating radar for determining volumetric soil water content;
111
result of comparative measurements at two sites. Journal of Hydrology, 316 – 338. Ward, S.H., 1990. Resistivity and induced polarization methods. Geotechnical and Environmental Geophysics, vol. 1. Society of Exploration Geophysicists, Tulsa Oklahoma, pp. 147 – 189. Yacoub, N.K., Scott, J.H., Mckeown, F.A., 1970. Computer ray tracing through complex geological models for ground motion studies. Geophysics 35 (4), 586 – 602. Yilmas, O., 1987. Seismic data processing. Society of Exploration Geophysicists, Tulsa Oklahoma. 526 pp.
Water content and porosity estimated from ground-penetrating radar and resistivity Anita Turesson Department of Geology, Earth Sciences Centre, Go¨teborg University, Box 460, SE-405 30 Go¨teborg, Sweden Received 1 September 2004; accepted 29 April 2005
Abstract Both ground-penetrating radar and the resistivity method have proven to be useful tools for exploring water content variations, since related parameters such as dielectric constant and the resistivity of rocks and sediments are highly dependent on moisture content. These methods were used independently to estimate volumetric water content in the unsaturated zone and porosity in the saturated zone in a 100-m sandy section. Two sample sites along the profile were also chosen for a shallow geophysical investigation and soil sampling, to enable the calibration and verification of the indirect geophysical methods. The grain distribution at these sites is dominated by medium-sized sand (0.25–0.5 mm). The water content was 6.9 vol.% and calculated porosities are 37% and 40% respectively. At each of these sites the mean water content values calculated from resistivity are within one percentage unit of measured water content while those calculated from ground-penetrating radar give higher values by as much as 2.9 percentage units. The water contents in the unsaturated zone in the section, estimated from resistivity and ground-penetrating radar, show very similar trends, although that deduced from ground-penetrating radar is generally somewhat larger, consistent with the results from the sample sites. The mean porosity values obtained from the two methods in the saturated zone are in good agreement. D 2005 Elsevier B.V. All rights reserved. Keywords: Ground-penetrating radar; Resistivity; Volumetric water content; Porosity
1. Introduction Soil water content and porosity are important variables in hydrological processes and are of primary interest in hydrogeological investigations. Ground-penetrating radar (GPR) has proven to be a
E-mail address: [email protected]. 0926-9851/$ - see front matter D 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.jappgeo.2005.04.004
promising technique for estimating water content in soil (Greaves et al., 1996; Van Overmeeren et al., 1997; Huisman et al., 2001). Using GPR in combination with mixing formulae, the water content can be estimated from dielectric constants, which are calculated from interval velocities of radar waves. The GPR method has successfully been applied to shallow (less than 50 m) geological surveys (Davis and Annan, 1989). The advantage of the method is
100
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
its vertical and lateral resolution in high resistivity environments. The main restriction of the method is the limited penetration in conductive materials, such as clays or soils with saline or contaminated pore water. Empirical relationships relating electrical resistivity to porosity have long been known (Archie, 1942). The water content in a clean sand formation can be estimated with Archie’s saturation equation. The applicability of geoelectric methods is not restricted to certain sedimentary sequences. However, the principle of equivalence, i.e. the thickness and resistivity of a unit can vary within certain limits and still give equivalent models (Parasnis, 1997), makes it difficult to estimate true subsurface resistivities. GPR and geoelectric methods have been used in combination to estimate water content and porosity in different ways. Dannowski and Yaramanci (1999) used GPR to constrain the geometry used in geoelectric inversion and compared the results from the two methods. Garambois et al. (2002) combined GPR and geoelectrics to estimate water content and water conductivity variations in the unsaturated zone. In contrast to previous investigations, this study aims to independently evaluate two methods, GPR and resistivity, used to assess water content (porosity in the saturated zone) according to the relationship of Topp et al. (1980) and Archie’s saturation equation (Ward, 1990) respectively. The Topp et al. (1980) equation, relating the dielectric constant to water content, was chosen among other mixing formulae (Mavko et al., 1999) because of its simplicity and good results (Greaves et al., 1996). Using standard techniques, in combination with well known empirical relationships, additional parameters can easily be obtained. This study uses geophysical data collected in a 100 m profile during two consecutive days in stable weather conditions. Two spots along the profile were chosen for small-scale shallow investigations to allow calibration and verification of the indirect geophysical methods. These investigations included GPR measurements, vertical electrical soundings, and soil sampling at a depth of 1 m. The measurements and sampling were made during two consecutive days at each site within a month of the profile survey. The data from each site were compared separately. The soil samples were analyzed to determine grain size distribution and gravimetric
water content, which was converted to volumetric water content. In addition, the P-wave seismic refraction was used to discriminate between the unsaturated and saturated zone, to support the results obtained from the GPR and resistivity. The objective of this paper is to compare and evaluate the variations in water content (porosity in the saturated zone) in sand estimated independently by two standard geophysical techniques, ground-penetrating radar and resistivity, in combination with the relationship of Topp et al. (1980) and Archie’s saturation equation (Ward, 1990).
2. Geological setting The test site is located at Veddige 70 km south of Go¨teborg on the Swedish west coast at an elevation of 15 m (Fig. 1). The post-glacial marine limit in the region is 65 m above sea level (Pa˚sse, 1986). During the overall post-glacial regression a small transgression (the Tapes transgression) occurred in certain parts of southwestern Sweden, reaching 17 m above sea level in the area (Pa˚sse, 1986). Large quantities of mostly glaciofluvial sediments, originating from a terminal moraine in the vicinity, were redeposited in an old channel. The area is underlain by wave-sorted sand and gravel which overlies the
Fig. 1. Map of southwestern Sweden showing the location of the site discussed in the text.
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111 Table 1 Grain size distributions presented as the mean of two samples at each site Sample site
Grain size distribution (mm)
Weightpercent
Measured h
Calculated porosity
30 m
Gravel (N2) Sand (0.071–2) Silt (b0.071) Gravel (N2) Sand (0.071–2) Silt (b0.071)
16 80 4 21 75 4
6.9%
39.9%
6.9%
37.3%
60 m
101
m. The information from these sample pits is summarized in Table 1.
3. Seismic refraction 3.1. Method
Measured volumetric water content (h) is the mean of three samples. A grain density of 2.65 g/cm3 was used in the porosity calculation. The samples were taken at about 1 m depth. See Figs. 3 and 7 for locations of sample sites.
glacial marine clay that in turn overlies sand and gravel (Pa˚sse, 1986). The surveyed profile is located on the redeposited sand with a maximum depth of 14 m. Soil samples were taken at a depth of 1 m at 30 and 60 m along the profile. These soil samples were analyzed for gravimetric water content and converted to volumetric water content using soil and pore water density estimates between 1.58– 1.68 and 1.00 g/cm3 respectively. The two sample sites gave the same result of 6.9 vol.%. The calculated porosities, using a grain density of 2.65 g/ cm3, are 37% and 40%. The pits revealed a sharp boundary between soil and sand at 0.33 m. The sand was analyzed to determine particle size distribution and is mostly of medium size (0.25–0.5 mm), 58% at sample site 30 m and 42% at 60
To determine the depth to the water table by a third independent method, we made a P-wave refraction survey using an ABEM Terraloc Mark 6 seismograph. The data were collected using 36 geophones (10 Hz) with 2-m spacing. Two spreads were measured, making a total length of 142 m. The energy source was a sledgehammer hitting a steel plate and three blows with the hammer were stacked in each record. The data were processed using the delay-time method (Pakiser and Black, 1957) followed by ray-tracing (Yacoub et al., 1970). 3.2. Result The seismic model shows the interface between unsaturated and saturated zones declining from 8.1 to 10.1 m in the investigated area (outlined in Fig. 2). The first layer with 375 m/s velocity is dry sand. The second layer with 1465 m/s velocity is interpreted as wet sand. The P-wave velocity increases greatly when water saturation reaches 100% (Bachrach and Nur, 1998), the refraction method is therefore suitable to determine the water table depth in coarser material such as sand. The water table deduced from seismic
Fig. 2. The seismic refraction model showing the interface between unsaturated and saturated zones. The water table declines from 8.1 to 10.1 m in the outlined area which coincides with the GPR profile.
102
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
refraction is plotted in all further models calculated from GPR and CVES.
4. Ground-penetrating radar 4.1. Method Two different GPR techniques were used in this study, the common-offset and the common midpoint (CMP) methods. The former was used to get an overview of the subsurface and to reveal any steeply declining surface, while the latter was used to get velocities with depths at regular intervals along the profile. Additional CMP soundings were measured at the two sample sites at 30 and 60 m to enable comparisons between water contents deduced from GPR and known water contents. GPR data is usually collected with the commonoffset method, in which the transmitting and receiving antennae are a fixed distance apart (Davis and Annan, 1989). The 100 m common-offset profile was measured using a Ramac GPR system from Mala˚ Geoscience. Antennae were 2 m apart with a centrefrequency of 50 MHz. The antennae were orientated parallel to each other and perpendicular to the profile. The data were collected every 0.2 m at a sample frequency of 500 MHz. To improve the signal-tonoise ratio, every trace was vertically stacked using the auto-stack option, which means that every trace was stacked between 32 and 64 times. A bhip chainQ, calibrated to actual length, was used to measure distance. This is a simple way to control the distance between each trace measured. A cotton-thread runs out as the operator slowly walks the profile. Every 0.2 m (in this case) the radar is triggered and a radar pulse is transmitted and subsequently received by the antenna. A more labour-intensive way to obtain data is to make an entire GPR survey using the multi-offset CMP technique (Fisher et al., 1992). In this case the radar data can be used for traditional seismic processing and are sorted in CMP-gathers, which are used in velocity analyses. The one-dimensional velocity models analyzed from CMP-gathers are then interpolated for the construction of a two-dimensional velocity profile. However, a less laborious way, used in this survey, to obtain such a 2-D velocity profile is to conduct individual multi-offset CMP soundings at
appropriate intervals along the profile. These measurements can easily be collected in the field using the hip chain measuring device. In a multi-offset CMP sounding the separation between transmitting and receiving antenna is continually increased from a fixed central location while the two-way travel time to reflectors are measured. Any subsurface contrast in electromagnetic properties results in energy being reflected back to the surface. Each reflection measured in this manner is used to derive the RMS (root mean square) velocity down to it. The CMP soundings were conducted using 50 MHz antennae with a sampling frequency of 500 MHz. At each CMP location the antenna separation was increased from 0 to 20 m, with increments of 0.2 m. The measuring device was placed at the midpoint. The true distance walked was corrected in the processing of data. To ensure that the antennae were moved equally from the midpoint, the measurements were performed stepwise using a measuring tape to control distances. Two persons are needed for this procedure. The spacing of individual CMP soundings, which should be measured in a profile, is a compromise between lateral and vertical variations of radar wave velocity and the time and effort to make the measurements. In this case we chose to conduct a spacing of 5 m and a total number of 21 CMP’s, which took about 2 h to collect. At the sample sites, CMP soundings were conducted in two directions to reveal possible three-dimensional geometry. The antennae centre frequency used was 200 MHz with a sampling frequency of 2000 MHz. The trace increment was 0.1 m. The vertical auto-stack function was used in all CMP measurements. The processing of GPR data (50 MHz) included time-zero adjustments and low-cut filtering (dewow), which removes low-frequency induction effects on the radar equipment. The data were also compensated for geometrical spreading and attenuation. The linear part of the gain was set to 0.02 (1/pulse width) and the exponential part was set to 0.03 dB/m (Davies and Annan, 1989). AGC scaling was used for display. The semblance approach (Yilmas, 1987) was used to pick preliminary RMS (or normal move-out) velocities. If the CMP data contains many and closely-spaced (in time) reflections it could be difficult to distinguish between real reflections arising from the interface between two electrically different media, and just a
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
complicated reflected wavelet. The compensation for geometrical spreading and attenuation helped to pick strong reflections in true amplitude display. To refine the velocity picks, hyperbolae were superimposed on the actual CMP gather to attain optimal fit. 4.2. Water content deduced from GPR From CMP-gathers the RMS velocity to reflectors is determined. The interval velocity, between reflectors, is calculated using the Dix (1955) equation: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v2RMS;n tn v2RMS;n1 tn1 Vi ¼ ð1Þ tn tn1 Where Vi is the interval velocity, v RMS are the RMS velocity and t n and t n 1 is the reflected ray two-way travel times to the nth and (n1)th reflectors respectively. The GPR produces high frequency electromagnetic energy in the 10 to 1000 MHz range. The term used to describe the high frequency electromagnetic properties of materials is called the dielectric constant or the relative permittivity. The complex dielectric constant is given by: KT ¼ KV þ jfKW þ ðrdc =xe0 Þg
ð2Þ
Where KV is the real part of the dielectric constant, and KW is the imaginary part of the dielectric constant or the electric loss, r dc is the zero-frequency conductivity, x is the angular frequency, e 0 is the free-space permittivity (8.854 10 12 F/m), and j is ( 1)1 / 2. In soils where the electric loss is small, that is with a lower conductivity than 10 mS/m, the following
103
relationship can be used to calculate the real part of the dielectric constant (KV; Davis and Annan, 1989): c ð3Þ V c pffiffiffiffiffi KV Where V is the propagation velocity of an electromagnetic wave in a medium with a real dielectric constant of KV, c is the velocity of an electromagnetic wave in free space (3 108 m/s). Topp et al. (1980) found that the real part of the dielectric constant (KV) seems to be highly sensitive to volumetric water content, but only weakly sensitive to soil type and density. They used a wide range of soil samples, sandy loam to clay, to derive an empirical relationship between the apparent (measured) dielectric constant and volumetric water content: hv ¼ 5:3 102 þ 2:92 102 Ka 5:5 104 Ka2 þ 4:3 106 Ka3
ð4Þ
Where h v is the volumetric water content (the ratio of water volume to total sample volume). For low-loss materials K a c KV where K a is the apparent dielectric constant. The water content (h) equals the product of porosity (/) and water saturation (S w). In water saturated soils the water content (h) is a measure of porosity (/). h ¼ /dSw
ð5Þ
4.3. Results The 100 m long, GPR common-offset profile is presented in a time–depth section in Fig. 3. The
Fig. 3. The GPR common-offset profile with depth represented by time. The arrows show the water table between 20 and 100 m at about 150 ns and the lowest picked reflector used in velocity analyses between 215 and 260 ns. Sample locations are indicated.
104
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
Fig. 4. Examples of CMP measurements used in the velocity analyses with a one-dimensional velocity model (continuous line: interval velocity, dotted line: RMS velocity), hyperbolic adaptions and semblance images at a) 30 m, b) 60 m and c) 90 m along the profile.
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
105
Fig. 4 (continued).
distinct reflector between 20–100 m at about 150 ns is interpreted as the water table in sand, which also is supported by the seismic refraction model (Fig. 2). The slightly declining reflector below that, from 215 to 260 ns, marks the lowest reflector picked in velocity analyses. Note the diffractions from cables (known
locations) in the upper part at about 72 and 85 m. Examples of CMP gathers with semblance and calculated interval velocities are shown in Fig. 4. The 21 CMP gathers were used to construct the two-dimensional velocity section converted to depths (Fig. 5). This section shows interval velocities between 72 and
Fig. 5. The GPR two-dimensional interval-velocity profile constructed from 21 one-dimensional velocity models. The continuous line shows the water table deduced from seismic refraction.
106
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
Fig. 6. Volumetric water content section calculated from GPR interval velocities using the relationship of Topp et al. (1980). The continuous line shows the water level deduced from seismic refraction.
149 m/As. The line marks the water table deduced from the seismic method. The mean interval velocity in the saturated zone is 76 m/As. The dielectric constants were calculated using Eq. (3) and the Topp et al. (1980) relationship (Eq. (4)) was used to convert the dielectric constants to volumetric water content (Fig. 6). The main part of the unsaturated zone shows water content between 10% and 12% but there are areas with lower water content (7–10%). The western part deviates with higher water content. In the saturated zone the water content varies between 22% and 31% (mean 28%). The CMP soundings at the sample sites revealed no three-dimensional heterogeneity. The interval velocities 127 m/As (at 30 m) and 125 m/As (at 60 m) yielded water contents of 9.4 and 9.8 vol.% respectively. These water contents differ by up to 2.9 percentage units from measured water content at these sites (Table 2).
5. Resistivity 5.1. Method The resistivity method is based on measuring the electrical potential which results from an applied direct electrical current flowing in the ground. The distribution of the electrical potential field depends in turn on the resistivity of the ground. In a multi-electrode array the measured result is displayed as a two-dimensional variation of apparent resistivity. Software applications for inversion of two-dimensional apparent resistivity to solve for true resistivity can be classified as either smooth inversion (DeGroot-Hedlin and Constable, 1990) or block inversion (Inman, 1975) methods, each of which has some disadvantages. Smooth inversion has a tendency to smear both resistivity and depth to interfaces even in the case of well-defined structures with sharp resistivity contrasts. On the other hand
Table 2 GPR wave velocity and vertical electrical soundings (VES) measured north–south (N–S) and east–west (E–W) directions Sample site
GPR velocity m/As
Calculated h (%)
VES q (Vm)
RMS error (%)
Calculated h (%)
30 30 60 60
127 127 125 125
9.4 9.4 9.8 9.8
11,300–11,500 –11,700 7500–7600–7800 6500–6600–6800 5600–5800–6100
0.9 1.3 1.6 3.2
5.3–5.3–5.2 6.8–6.7–6.6 7.4–7.3–7.2 8.1–7.9–7.7
m (N–S) m (E–W) m (N–S) m (E–W)
Equivalent VES models resistivities (q) up to 1.2% fit, and the RMS (root mean square) errors for the best model (middle value) are given. To calculate the volumetric water content (h) from Archie’s formula (Eq. (7)), pore water conductivity of 13 mS/m and m = n = 1.7 was used. See Figs. 3 and 7 for locations of sample sites.
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
block inversion requires a starting model close to the truth, which is rarely known precisely and is difficult to construct especially in complex cases. We decided to use the smooth inversion routine. The continuous vertical electrical soundings (CVES) were conducted with an ABEM Lund imaging system, using Wenner measurement geometry. The electrode spacing varied from 3 to 72 m in the 240 m profile. The CVES profile was centred on the GPR profile. A software package based on the smoothnessconstrained least-square method (DeGroot-Hedlin and Constable, 1990) was used to invert apparent resistivity to true resistivity. At the sample sites, 30 and 60 m, vertical electrical soundings (Schlumberger geometry) were measured in two directions, using the ABEM equipment. A total of 13 measurements were done moving the current electrodes equally and stepwise from the midpoint, with distances increasing from 0.5 to 8 m between the midpoint and current electrodes. 5.2. Water content deduced from resistivity Electrical conduction in soil is largely electrolytic, taking place in connected pore spaces and along grain boundaries. The relationship between resistivity (the inverse of conductivity) and porosity in sedimentary clay-free rocks is expressed by the formation factor, which is the ratio of the resistivity of the porous media to that of the pore fluid (Archie, 1942; Ward, 1990). q F¼ ¼ ad/m ð6Þ qw Where F is the formation factor, q is the bulk resistivity of the rock, q w is the resistivity of the
107
pore fluid, / is the porosity, and a and m are constants. A general form of Archie’s saturation equation is: q ¼ qw /m Swn
ð7Þ
Where q and q w is the bulk resistivity of the rock and the resistivity of the water respectively, / is the porosity, S w is the fractional water saturation, and n is the saturation exponent, which normally is equal to 2 (Ward, 1990). For a water-saturated rock Eq. (7) is reduced to: q ¼ qw /m
ð8Þ
Jackson et al. (1978) found that the exponent m was dependent on the shape of the particles, increasing as they became less spherical, while variation in size appeared to have little effect. Samples of natural sand have values of m in the range 1.4 to 1.6. In solving Archie’s formula (Eq. (7)) for volumetric water content, which is the product of porosity and water saturation, parameters such as the electrical conductivity of pore water, m and n had to be estimated. To do this, volumetric water content was calculated using a range of conductivities (2–40 mS/ m) and m = n(1.3–2). After comparison between calculated and measured water content from the two sample sites conductivity was estimated at 13 mS/m and m = n = 1.7. This conductivity for pore water is supported by measurement in a nearby well which was also 13 mS/m. In this investigation it was assumed that these parameters would not change in the section.
Fig. 7. The resistivity model. The outlined area coincides with the GPR profile and is used in the water content calculation. Sample locations are indicated.
108
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
5.3. Results
6. Discussion and conclusion
The 240 m long resistivity profile is presented in Fig. 7. The area investigated and compared to GPR measurements is outlined. Two anomalies with lower resistivities are seen at about 72 and 85 m. These are probably due to artefacts from cables, also seen in the GPR section (Fig. 3). Although the cables are less than 1 m below surface, the effects of resistivity distortion are seen to considerably greater depths. The water content was calculated using Archie’s equation (Eq. (7)), with water resistivity set to 77 Vm (the inverse of 13 mS/m) and m = n = 1.7. The resistivitybased water content section (h resistivity; Fig. 8) is more varied than the corresponding GPR section (h GPR; Fig. 6), containing both lower and higher values, although the general trend is the same with the exception of the two artefacts mentioned above. The water content in the unsaturated zone is between 5% and 14% with the lowest value about 2 percentage units lower than in corresponding h GPR profile. The increasing water content in the western part can also be seen in the h GPR profile. In the saturated zone the water content varies between 16% and 40% (neglecting the anomalies). The electrical soundings at the sample sites show some anisotropy, especially at 30 m. At each site mean water content values calculated from resistivity, 6.0% (30 m) and 7.6% (60 m), are within 1 percentage unit of the measured water content, 6.9% (Table 2).
Referring to data from the sample pits (Tables 1 and 2) the volumetric water content calculated from the Topp et al. (1980) equation is higher than both the measured value and that calculated from the Archie equation. Considering the corresponding profiles, h GPR and h resistivity (Figs. 6 and 8), the same trend is true for the larger part of the section, which can be seen in Fig. 9 showing the difference between these two results. However, except for the two anomalies the h resistivity shows higher water content near to the surface. This could be due to the fact that the smallest electrode spacing is 3 m which enables a shallower measurement of water content than for h GPR where the first picked reflector in the velocity analysis is at a depth of about 4 m. The h resistivity profile seems to be more detailed due to a denser sample grid with depth, but also probably due to the smooth inversion routine, which renders a gradual change of resistivities even if sharp boundaries exist in the subsurface (Olayinka and Yaramanci, 2000, 2002). On the other hand the less detailed h GPR is restricted by the number of reflections present in the subsurface and by the distance between the CMP’s. So in general, the difference in data density accounts for some of the differences between the two methods seen in Fig. 9. In the saturated zone the water content, which is the porosity when the pores are saturated, varies
Fig. 8. The volumetric water content section calculated from resistivity using Archie’s (1942) equation. Water conductivity is set to 13 mS/m and m = n = 1.7. The continuous line shows the water table deduced from seismic refraction.
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
109
Fig. 9. This section shows the difference between volumetric water contents deduced from GPR and resistivity.
between 22% and 31% in h GPR (Fig. 6) with a mean of 28%. The corresponding range for h resistivity (Fig. 8) is wider, roughly between 16% and 40% but also with a mean of 28%, calculated between depths of 10.80 and 12.89 m (excluding the anomaly between 82 and 90 m). In the sample pits the porosity was measured at 37% and 40% respectively. This suggests a decrease in porosity by approximately 10 percentage units at about 10 m depth. However, the measured porosity from 1 m depth is the total porosity which is the sum of effective, trapped and isolated porosities whereas the porosity calculated at 10 m depth from the resistivity method is the effective porosity because the current is largely electrolytic. The GPR method is based on wave propagation by analogy with seismic methods, so in that sense the porosity calculated from GPR is the total porosity. As the methods give very similar results this would indicate that the effective porosity equals the total porosity in this case and that no trapped and isolated pores exist. Porosity is governed by many factors such as the uniformity of grain size, (sorting), grain shape, packing, and compaction during and after deposition. Packing alone can contribute significantly to the difference in porosity. The end members of packing modes for spheres of uniform size, the cubic and rhombohedral packing, have porosities of 48% and 26% respectively (Graton and Fraser, 1935). Sorting also has a large influence on the porosity: up to 25% difference between well-sorted and very poorlysorted sands of the same mean grain size were reported by Beard and Weyl (1973). The dominant
factors in this study are not known but if the results indicated by two methods are correct, it is likely that more than one factor are responsible for the relatively large decrease in porosity with depth. Unfortunately, no direct control is available to verify the porosity at 10-m depth. A precondition for using the GPR in this type of study is the presence of several electromagnetically contrasting layers. This is not often a problem as sandy material is commonly stratified. Another requirement is that the layers are more or less horizontal since the equation of Dix (1955) is valid only for horizontal surfaces. The Topp et al. (1980) equation requires values for dielectric constant for each layer, derived from the interval velocity, as input for water content calculation. Thus the RMS velocity has to be carefully picked for each layer. In this study these were determined using both the semblance approach to select strong reflections and hyperbolae fitted to CMP gathers to refine the two-way times and RMS velocities. These velocities were thoroughly analysed and the maximum estimated error could be F1 ns in two-way-time and F0.003 m/ns in RMS velocity, which would result in about 1.5% error in water content. The smooth inversion method of DeGroot-Hedlin and Constable (1990), used for the geoelectrical data, has the advantage of being fully automatic and not needing any prior information, however the gradual change from high to low resistivities which is inherent to this method makes it difficult to determine an intrinsic value for porosity. A prerequisite for a reliable water-content model deduced
110
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111
from resistivity is that true resistivities can be derived. This could also be problematic due to the equivalence principle relating depth and resistivity. An alternative inversion routine, the block inversion (Inman, 1975), has proven to be useful when the subsurface layer geometry is simple (Dannowski and Yaramanci, 1999). In more complex environments the input model required for this inversion routine may be too difficult to estimate to get a satisfying result. In using the Archie equation several parameters, such as the constants m and n and sometimes also water conductivity, have to be estimated. The small-scale investigation at the sample sites was used to calibrate these parameters. Furthermore, it was assumed that these parameters did not change in the section. In this survey this assumption can be justified by the facts that the constants m and n are related to pore shape and pore fluid and that the wave-sorted sandy material in the section has the same origin and genesis as sand washed out from a terminal moraine. In summary, two methods, ground-penetrating radar and resistivity, were independently evaluated for their capability to assess water content and porosity in a sandy section. The methods were used in combination with empirical relationships. Additional information from two sample sites was used to relate the indirect methods to known water content and porosity. This is of special importance when using empirical relationship in different specific environments. The results obtained showed very similar trends of water-content distribution, although absolute values differ somewhat, and there is a good agreement between the methods in the saturated zone if the mean porosity is compared. The use of two independent methods greatly strengthens the results which can be obtained in this type of study.
Acknowledgement I would like to thank associate professor Gustaf Lind, Earth Sciences Centre, Go¨teborg University, for assistance during fieldwork and for constructive review of the manuscript. I also thank professor David Cornell, Earth Sciences Centre, Go¨teborg University, for correction of the English language.
References Archie, G.E., 1942. The electrical resistivity log as an aid in determining some reservoir characteristics. Transactions of the American Institute of Mining and Metallurgical Engineers/Petroleum Division 146, 54 – 62. Bachrach, R., Nur, A., 1998. High-resolution shallow-seismic experiments in sand: Part I. Water table, fluid flow, and saturation. Geophysics 63 (4), 1225 – 1233. Beard, D.C., Weyl, P.K., 1973. The influence of texture on porosity and permeability on unconsolidated sand. The American Association of Petroleum Geologists Bulletin 57, 349 – 369. Dannowski, G., Yaramanci, U., 1999. Estimation of water content and porosity using combined radar and geoelectrical measurements. European Journal of Environmental and Engineeering Geophysics 4, 71 – 85. Davis, J.L., Annan, S.P., 1989. Ground-penetrating radar for highresolution mapping of soil and rock stratigraphy. Geophysical Prospecting 37 (5), 531 – 551. DeGroot-Hedlin, C., Constable, S., 1990. Occam’s inversion to generate smooth, two-dimensional models from magnetotelluric data. Geophysics 55 (12), 1613 – 1624. Dix, C.H., 1955. Seismic velocities from surface measurements. Geophysics 20 (1), 68 – 86. Fisher, E., McMechan, G.A., Annan, A.P., 1992. Acquisition and processing of wide-aperture ground-penetrating radar data. Geophysics 57 (3), 495 – 504. Garambois, S., Se´ne´chal, P., Perroud, H., 2002. On the use of combined geophysical methods to assess water content and water conductivity of near-surface formations. Journal of Hydrology 259, 32 – 48. Graton, L.C., Fraser, H.J., 1935. Systematic packing of spheres, with particular reference to porosity and permeability. Journal of Geology 43, 785 – 909. Greaves, J.G., Lesmes, D.P., Lee, M.L., Tokso¨z, M.N., 1996. Velocity variations and water content estimated from multi-offset, ground-penetrating radar. Geophysics 61 (3), 683 – 695. Huisman, J.A., Sperl, C., Bouten, W., Verstraten, J.M., 2001. Soil water content measurements at different scales: accuracy of time domain reflectometry and ground-penetrating radar. Journal of Hydrology 245, 48 – 58. Inman, J.R., 1975. Resistivity inversion with ridge regression. Geophysics 40 (5), 798 – 817. Jackson, P.D., Taylor, S.D., Stanford, P.N., 1978. Resistivity–porosity particle shape relationships for marine sands. Geophysics 43 (6), 1250 – 1268. Mavko, G., Mukerji, T., Dvorkin, J., 1999. The rock physics handbook. Tools for Seismic Analysis in Porous Media. Cambridge University Press. 329 pp. Olayinka, A.I., Yaramanci, U., 2000. Use of block inversion in the 2-D interpretation of apparent resistivity data and its comparison with the smooth inversion. Journal of Applied Geophysics 45, 63 – 81. Olayinka, A.I., Yaramanci, U., 2002. Smooth and sharp-boundary inversion of two-dimensional pseudosection data in presence of a decrease in resistivity with depth. European Journal of Environmental and Engineering Geophysics 7, 139 – 165.
A. Turesson / Journal of Applied Geophysics 58 (2006) 99–111 Pakiser, L.C., Black, R.A., 1957. Exploring for ancient channels with the refraction seismograph. Geophysics 22 (1), 32 – 47. Parasnis, D.S., 1997. Principles of Applied Geophysics. Chapman & Hall, UK. 429 pp. Pa˚sse, T., 1986. Beskrivning av jordartskartan Kungsbacka SO. Sveriges Geologiska Underso¨kning Ae 56 (106 pp.). Topp, G.C., Davis, J.L., Annan, A.P., 1980. Electromagnetic determination of soil water content: measurements in coaxial transmission lines. Water Resources Research 16 (3), 574 – 582. Van Overmeeren, R.A., Sariowan, S.V., Gehrels, J.C., 1997. Ground penetrating radar for determining volumetric soil water content;
111
result of comparative measurements at two sites. Journal of Hydrology, 316 – 338. Ward, S.H., 1990. Resistivity and induced polarization methods. Geotechnical and Environmental Geophysics, vol. 1. Society of Exploration Geophysicists, Tulsa Oklahoma, pp. 147 – 189. Yacoub, N.K., Scott, J.H., Mckeown, F.A., 1970. Computer ray tracing through complex geological models for ground motion studies. Geophysics 35 (4), 586 – 602. Yilmas, O., 1987. Seismic data processing. Society of Exploration Geophysicists, Tulsa Oklahoma. 526 pp.