VZJ sign up for etocs
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (10)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Cassiani, G.
Right arrow Articles by Gallotti, L.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Cassiani, G.
Right arrow Articles by Gallotti, L.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Cassiani, G.
Right arrow Articles by Gallotti, L.
Related Collections
Right arrow Water Content
Right arrow Ground Penetrating Radar, GPR
Right arrow Inverse Procedures/Parameter Estimation
Published in Vadose Zone Journal 3:1093-1105 (2004)
© 2004 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

SPECIAL SECTION: HYDROGEOPHYSICS

Vertical Radar Profiles for the Characterization of Deep Vadose Zones

Giorgio Cassiania,*, Claudio Strobbiab and Laura Gallottic

a Dipartimento di Scienze Geologiche e Geotecnologie, Università di Milano Bicocca, Italy
b European Centre for Training and Research in Earthquake Engineering, Pavia, Italy
c Dipartimento di Scienze dell'Ambiente e del Territorio, Università di Milano Bicocca, Italy
* Corresponding author (giorgio.cassiani{at}unimib.it)

Received 1 February 2004.



   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 MATERIALS AND METHODS
 MODELING
 RESULTS
 CONCLUSIONS
 REFERENCES
 
The deep vadose zone, down to a few tens of meters below the soil surface, is difficult to investigate and characterize, especially from a hydraulic point of view. We present a method well suited for this purpose, that makes use of a radar transmitting antenna at the surface and a radar receiving antenna deployed in a plastic-cased borehole. The technique is often referred to as vertical radar profiling (VRP). Vertical radar profiling yields information on (i) the distribution of radar velocity as a function of depth, from which moisture content distributions can be inferred, and (ii) the presence of reflecting horizons within the formation, often associated with lithological contacts. High-resolution, time-lapse VRPs were acquired from October 2002 to December 2003 at a contaminated site in Trecate, Northern Italy, from several existing boreholes originally installed for monitoring and remediation, with the aim of characterizing the dynamic hydrologic behavior of the deep vadose zone. We discuss data acquisition, processing, and inversion of these VRP data to derive moisture content profiles as a function of time. The VRP-derived moisture content profiles were used to calibrate a dynamic Richards' equation model via a Monte Carlo inversion approach. In this way, it was possible to identify the value of subsurface hydraulic parameters, in particular hydraulic conductivity, of the main lithological units, and the parameters' uncertainty. The location of lithology contacts was derived via analysis of up-going radar reflection events contained in the same VRP data.

Abbreviations: b.g.l., below ground level • GPR, ground penetrating radar • VRP, vertical radar profiling


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
MATERIALS AND METHODS
 MODELING
 RESULTS
 CONCLUSIONS
 REFERENCES
 
THE IDENTIFICATION of flow and transport characteristics of the vadose zone is a fundamental step in the characterization of contaminated sites and the estimation of the resulting risk of groundwater pollution. Not many techniques are capable of providing information regarding the conditions of deep vadose zones (deeper than several meters below ground level [b.g.l.]). Cross-borehole radar has gained popularity in this context for monitoring moisture content changes, thanks to its apparent simplicity and its high resolution capabilities (Lager and Lytle, 1977; Daily and Lytle, 1983; Eppstein and Dougherty, 1998; Hubbard et al., 1997; Parkin et al., 2000; Trinks et al., 2001; Binley et al., 2001, 2002a; Alumbaugh et al., 2002; Day-Lewis et al., 2002, 2003; Ferré et al., 2003; Tronicke et al., 2004). Cross-borehole radar tomography consists of sending a sequence of high-frequency electromagnetic signals through a transmitting antenna and measuring signals at a receiving antenna. One or both antennas are located in a plastic-cased borehole. The depth of antennas is progressively changed, thus exploring the velocity, and possibly the attenuation, characteristics of radar signal propagation along different directions in the subsurface.

The most common application of radar for hydrologic purposes is the identification of dielectric properties of the medium via measurement of travel time of the signal from transmitter to receiver. In low loss materials and at high frequencies (10–1000 MHz), the electromagnetic signal velocity depends on the dielectric properties only. The relative permittivity, or bulk dielectric constant {kappa}r, is derived from

[1]
where c is the radar wave velocity in air ({cong}0.3 m ns–1), and v is the measured radar velocity in the medium. The dielectric properties of a partially saturated porous medium are strongly linked to its moisture content because of the high dielectric constant of water ({kappa}r {cong} 81). The link between moisture content and dielectric constant has been the subject of research for more than two decades (Topp et al., 1980; Wharton et al., 1980; Roth et al., 1990; Knoll and Knight, 1994; Pham, 2000) and several semiempirical and empirical models are available for the conversion. Even though some uncertainty remains, these relationships are now widely used for in situ moisture content measurement in shallow soil, using techniques such as time domain reflectometry (TDR) (Topp and Davis, 1985; Topp et al., 2003) and frequency domain reflectometry (Dirksen and Hilhorst, 1994). The use of radar, particularly in cross-hole mode, for hydrologic purposes is based on the same physical principles and has become increasingly popular (for a review, see Huisman et al., 2003).

The obvious advantages of radar over other methods (such as TDR or neutron probe) include the following:

The use of radar in boreholes has tremendous advantages over surface to surface radar applications, especially since the depth of penetration is not drastically limited by conductive layers close to the surface, as in the case of surface deployed antennas, and not as much resolution is lost with depth as is observed for surface applications. An obvious disadvantage of borehole applications is the need for (plastic cased) boreholes. A further disadvantage of hole to hole radar is that the boreholes must be closely spaced (a few meters to a few tens of meters, under ideal conditions), a feature rarely available at sites of practical interest. Unlike cross-hole applications, VRPs (e.g., Zhou and Sato, 2000; Knoll and Clement, 1999; Clement and Knoll, 2001; Buursink et al., 2002) require only one borehole, with practical and financial benefits. The transmitting antenna is then placed on the ground surface while the receiving antenna is lowered into the borehole. The receiver, rather than the transmitter, is lowered into the holes to reduce the effect of ambient electromagnetic noise during the acquisition. However, VRP has one obvious disadvantage with respect to hole to hole applications: the need to extract interval radar velocity values from integral information of travel time from surface to depth. In fact, the acquisition geometry of VRP requires that data be inverted in a more sophisticated manner than commonly needed for zero-offset cross-hole radar profiles (e.g., Binley et al., 2001). In this respect VRP is similar to multiple offset gather tomography. Also, because with VRP one antenna is located at the ground surface, penetration may still be limited by the presence of shallow conductive layers, and vertical resolution is potentially affected by the change in signal frequency content as the receiver antenna is lowered deeper into the subsurface. However, both phenomena are much less dramatic than in surface to surface radar surveys. In particular, VRP resolution is more linked to accurate depth position of the receiver antenna than to the criteria of wavelet interference, such as Fresnel zone considerations, often used in surface radar and seismic applications.

The moisture content data derived from cross-hole radar and VRP can be used in inverse hydrogeophysical models. The aim is to characterize the hydraulic properties of the subsurface by constraining predictive models (e.g., based on Richards' equation) on available geophysical and hydrological data. It has been proven that the use of cross-hole radar information in time-lapse mode can improve such inversion dramatically (Hubbard et al., 2001; Binley and Beven, 2003; Binley et al., 2002b, 2004; Kowalsky et al., 2004), albeit the practical utilization of the approach is still in its infancy. No such exercise has been performed, to our knowledge, using ground penetrating radar (GPR) data in VRP configuration. Also, no attempt has been made so far to introduce VRP reflection information into such a conceptual framework.

Consequently, our objectives were:


   MATERIALS AND METHODS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
MODELING
 RESULTS
 CONCLUSIONS
 REFERENCES
 
Field Site
High-resolution, time-lapse VRPs were acquired at a crude oil–contaminated site in Trecate, Piemonte, Northern Italy, using a few existing boreholes originally developed for groundwater monitoring and remediation via bioventing (Fig. 1) . In 1994 the site was the scene of an inland crude oil spill following an oil well blowout (Well Tr 24). Details of the incident, which resulted in approximately 15000 m3 of middleweight crude oil being released at the soil surface and contaminating soil, vadose zone and groundwater, together with descriptions of subsequent site remediation were reported elsewhere (Reisinger et al., 1996; Brandt et al., 2002; Christensen et al., 2004). The Po river plain aquifer at Trecate comprises an extensive, unconfined silty sand and gravel unit more than 60 m thick beneath the site. The hydraulic properties of the aquifer in the Trecate region were determined from pumping tests performed in 1994 using six wells. The test data yielded average horizontal hydraulic conductivity, transmissivity, and porosity values of 56.5 ± 5.1 m d–1, 2700 ± 240 m2 d–1, and 0.29, respectively. The vertical hydraulic conductivity in the saturated zone was estimated to be 3.4 m d–1 (Burbery et al., 2004). Groundwater levels at the site show seasonal fluctuations of 5 to 6 m, with higher levels experienced during the summer as a result of surficial recharge from irrigation. Groundwater at the site flows in an easterly direction toward the Ticino River a few kilometers from the site. Information on site stratigraphy can be derived from many sources, including drilling logs, Geoprobe (Salina, KS) cores extracted for contamination monitoring across the site, and visual inspection of a nearby gravel quarry. However, no characterization via detailed geophysical well logs was ever attempted at the site. A thick sequence of poorly sorted silty sands and gravels in extensive lenses, typical of braided river sediments, characterizes the site stratigraphy. An artificial layer of clayey-silty material, <1 m thick, originally placed as a liner for rice (Oryza sativa L.) paddies overlies most of the site. Periodic VRP acquisition was performed in five different boreholes (Fig. 1), with approximately one survey every 2 wk for more than 1 yr, from October 2002 to December 2003. In this discussion we focus on the data collected from two boreholes: Borehole FW16 from December 2002 to February 2003 (during water table fall) and Borehole BB from March 2003 to November 2003 (during water table rise). The FW boreholes had to be removed in February 2003 as part of the site restoration plan. They were drilled using a 25.4-cm (10 inch) drill bit and equipped with 10.16-cm (4 inch) plastic casing and gravel backfill. Of the B boreholes, used for groundwater monitoring and drilled using a 10.16-cm (4 inch) bit and equipped with 5.08-cm (2 inch) plastic casing, only BB was permanently accessible during summer, when a crop was growing on the land. In spite of the different drilling size and completion characteristics, the above boreholes delivered data of excellent quality. However, the presence of the large gravel pack may have altered the hydrological meaning of data collected in the FW boreholes, since the VRP surveys may have been more sensitive to the moisture conditions in the gravel pack than in the surrounding formation. For this reason we concentrate mostly on the results obtained from borehole BB.



View larger version (58K):
[in this window]
[in a new window]
  		Fig. 1. Map of the Trecate site with vertical radar profiling (VRP)-monitored boreholes. Tr 24 is the well that caused the oil spill in 1994.

 
Methodology
Acquisition and Processing
The acquisitions were performed using a PulseEKKO 100 system (Sensors & Software Inc., Mississauga, ON, Canada) with 100-MHz borehole antennas. The transmitter antenna was placed radially with respect to the borehole axis, its center positioned at 0.70 m from the borehole axis. The receiver antenna was lowered in 5-cm increments, starting at 0.75 m below ground. For each point at depth, the signal was stacked 64 times to improve the signal/noise ratio. The adopted time sampling interval was 0.4 ns for all surveys. Periodically, a test was made on the VRP signal repeatability by placing the transmitting antenna on the soil surface radially along different azimuths, and checking the different arrival times at the receiver antenna. Note that this procedure does not alter the configuration of antennas in terms of mutual polarization. However, there is a potential problem in that anisotropy along different azimuth directions may be measured instead of signal repeatability. Given the nature of the site geology (i.e., sand-gravel quaternary sediments), there is no reason to believe that this type of anisotropy is present. For other geological conditions, repeatability could be assessed by repositioning the transmitter antenna in the same position and reacquiring the entire sequence of shots with the receiver antenna at depth. If azimuth rotation of the transmitter is allowed, as in our case, a further advantage is that the resulting error estimation incorporates implicitly also (minor) discrepancies between the site conceptualization (as a purely one-dimensional sequence of layers vs. depth) and reality, which might be more complex. We strongly suggest that some type of data quality assessment, such as the one proposed here, be made to evaluate the reliability and repeatability of the derived radar velocity estimates. An average error of roughly 0.5 ns was measured as long as the receiver antenna was located above the water table (as deep as 10 m b.g.l.), thus giving a good confidence on the quality of data collected in the vadose zone (Fig. 2) . When immersed below the water table, the signal repeatability decreases to about 4 ns. Therefore porosity estimates from below the water table (see Knoll and Clement [1999] and Buursink et al. [2002], who reported on sites with much shallower water table) are not reliable in our case. Note that a simple repetition of shots, as suggested by Lizarralde and Swift (1999) for seismic profiles, does not seem appropriate for VRPs. In seismics, one of the main sources of error is the difficulty of repeating shots with the same characteristics, be it with explosives, weight dropping, vibrators, or other sources. In the case of GPR measurements, repeatability errors arise mostly from antenna positioning (vs. depth and in the horizontal location in the borehole) and may result also from model errors (assuming the model is one-dimensional). Our choice of studying repeatability by considering complete independent surveys is consistent with existing literature on cross-hole radar repeatability. Indeed, the error estimates (around 0.5 ns) are very similar to the results of Alumbaugh et al. (2002), who also used a PulseEKKO borehole system with 100-MHz antennas. We also compared our error estimates with those computed using the discrepancy principle (Hansen, 1992) that is purely based on the characteristics of the regularized velocity inversion (see below), and found excellent agreement between the results of the two approaches. Note however that (i) application of the discrepancy principle requires some subjective choice of a "corner" on the computed error curve (for details see Moret et al., 2004), and (ii) such a principle is based entirely on the statistical properties of the regularized inversion, and therefore depends a priori on the type of (smooth) velocity model assumed to exist. Adopting an overall repeatability analysis is a more robust approach for data inversion. Note that a link exists between the travel time picking errors and the vertical resolution of VRP. As the receiver antenna is moved down the borehole in steps (e.g., 5 cm), the information is given by the travel time between transmitter and receiver. Therefore, the vertical spatial resolution depends mostly on our ability to distinguish travel time values of neighboring traces, which is strongly linked to the travel time picking error as discussed above. If travel times between two locations at depth differ less than the repeatability error for the adopted methodology (say, 0.5 ns), then the two depths are indistinguishable from the VRP viewpoint (i.e., the two locations fall into the same pixel in the "blurred" image of radar velocity given by VRP). To draw conclusions on the vertical resolution of our data, we computed the standard deviation of arrival times within vertical moving windows along the profile (Fig. 3) , and then compared this standard deviation with the absolute value of the travel time picking error measured in the repeatability analysis (Fig. 2). This approach shows that a reasonable estimate of vertical resolution for our VRP data is about 0.25 m, which corresponds to a standard deviation (above the water table) of about 0.5 ns. Note that this resolution is considerably smaller than the antenna length and the dominant signal wavelength (around 1 m), and larger than the vertical spatial sampling (5 cm). This scale of measurement is not dissimilar from the one derived by Buursink et al. (2002) (50 cm). Picking of first breaks was preferred to alternative approaches (such as picking of first peaks as in Buursink et al., 2002) because of the dispersive nature of the medium, especially at large depths. The wavelet changes shape because of higher frequency attenuation and propagation of different frequencies at slightly different velocities. This causes the wavelet to lengthen, while the time interval between first break and first peak changes with depth, thus causing a nonstationary error that is difficult to account for. It should be noted, however, that because of a decreasing signal/noise ratio with depth, when performing first break picking, there may be tendency to shift picks to later times. We deem that this error, in our case, is significantly less serious than those introduced by first peak picking. Vertical radar profiling data processing is relatively simple, involving:



View larger version (28K):
[in this window]
[in a new window]
  		Fig. 2. Example of error analysis on vertical radar profiling (VRP) data at the Trecate site. The black dots are the absolute values of arrival time differences obtained by rotating the transmitter antenna at the surface by 90° in azimuth. The gray continuous lines represent standard deviation of arrival times computed over moving windows of different size along the vertical direction. A comparison between moving window standard errors and difference in arrival times at different transmitter azimuth (around 0.5 ns) indicates that the actual resolution of VRP data at this site is about 25 cm. Data were collected at 5-cm vertical spacing.

 


View larger version (43K):
[in this window]
[in a new window]
  		Fig. 3. Computation of travel time error estimates using vertical moving windows: Borehole BB, 31 Mar. 2003. Note how the 1-m window, equal in size to the deployed antenna, overestimates the data error: within 1 m the trend in the data is clearly identifiable as a signal above the error level. On the contrary, arrival times within a window of 0.25 m are within the error range (0.5 ns).

 
Dewowing has proved critical to reliably estimate travel times from VRP records at Trecate. The problem caused by the "wow" becomes more severe as the receiver antenna is lowered deeper into the subsurface. This happens because, while the wow has the same energy at all depths, the radar wavelet signal becomes weaker with depth. Considering that the wow appears on the trace even before the first break of the radar wavelet, correct picking of the travel times at larger depth is problematic without dewowing. As pointed out by Gerlitz et al. (1993), the use of simple low-frequency removal filters may distort the wavelength shape, thereby introducing precursors before the first break. To preserve the timing of the first break, Gerlitz et al. (1993) suggested using a residual median filter, and provided relevant rules for calculating the median filter length. The choice of a suitable dewowing algorithm depends on how the (dewowed) data will be used. At least three possible applications can be envisioned: (i) velocity tomography, for which correct picking of first arrivals is necessary; (ii) reflection processing, for which preservation of the signal characteristics after the first break is essential; and (ii) attenuation processing, for which preservation of the spectrum characteristics after removal of the wow is necessary. Our experience shows that while dewowing via residual median filtering is a good approach for (i), it fails to provide good filtering characteristics for the other two applications, for which the traditional residual mean filtering approach is more suitable. Figure 4 shows an example of data from shallow and deep parts of one of the VRP surveys at Trecate, processed both via residual median and mean dewowing. According to the approach by Gerlitz et al. (1993), a filter length of 80 ns was adopted for the residual median filter. Figure 4 shows also results of applying residual mean filters of the same length (80 ns) and shorter length (10 ns), the latter being the best approach to preserve spectral characteristics (not shown). It is clear from Fig. 4 that residual median dewowing is the only suitable approach for first break picking. Figure 5 shows two examples of entire VRP profiles processed with either type of dewowing. Residual mean filtering produces better results in terms of preservation of reflected events, and has consequently been adopted for our VRP reflection analysis.



View larger version (31K):
[in this window]
[in a new window]
  		Fig. 4. Dewowing of vertical radar profiling (VRP) data from Trecate. The upper row of figures refers to data from 1.5 m b.g.l.; the lower row refers to 8.5 m b.g.l. Plots (a) and (e) are raw field data. (b) and (f) are the corresponding dewowed data using a residual median filter with length equal to 80 ns. (c) and (g) are dewowed data using a residual mean filter with length equal to 80 ns. (d) and (h) are dewowed data using a residual mean filter with length equal to 10 ns. The detrimental effect of the "wow" is more serious at depth, requiring that dewowing is applied (note that the "wow" in (e) makes first break picking impossible). The only dewowing algorithm that worked properly for all depths was the residual median filter (b) and (f). Spurious precursors were introduced by the residual mean filters.

 


View larger version (95K):
[in this window]
[in a new window]
  		Fig. 5. Examples of vertical radar profiling (VRP) data from the Trecate site. Time zero correction and trace normalization were applied to all datasets. Data in (a) and (b) were processed with residual median dewowing, 80-ns filter length, while the corresponding data in (c) and (d) were processed with residual mean dewowing, 10-ns filter length. The latter approach was more effective at preserving reflected events. Note the distinct slope change in the first arrivals in correspondence of the water table depth and the clear up-going reflections at several depths (particularly around 2 and 6 m b.g.l.). See Fig. 7 for the geometry of reflected events.

 


View larger version (30K):
[in this window]
[in a new window]
  		Fig. 7. Schematic representation of reflection mechanisms in vertical radar profiling (VRP) and the resulting zero-offset radargram.

 

   MODELING
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 MODELING
RESULTS
 CONCLUSIONS
 REFERENCES
 
Inversion
Travel times, as derived from picking of first arrivals, must be inverted into interval radar velocity values. The velocities can be transformed into estimates of moisture content via suitable petrophysical models (discussed in the section below). Such inversion involves as many unknowns as equations, since for each (say 5-cm) downward displacement of the receiving antenna, a new travel time is measured, and the corresponding interval velocity of the 5-cm-thick layer is the relevant unknown. This results in an even-determined problem:

[2]
where T is the vector of first-arrival times, P is the vector of (unknown) interval slownesses (reciprocal of velocities), and A is a lower triangular matrix defined as (Fig. 6)

[3]



View larger version (15K):
[in this window]
[in a new window]
  		Fig. 6. Scheme considered for vertical radar profiling (VRP) data inversion, from first arrival times to interval radar velocity values. In this study the thickness hi was 5 cm.

 
The simple formulation (Eq. [2]) is based on a straight ray approximation to the actual radar energy propagation (Fig. 6), valid as long as the ratio of source offset/station depth is <1 (Schuster et al., 1988). However, the apparent simplicity of Eq. [2] hides the inherent instability of the solution caused by errors in the data (travel times). If such errors are not accounted for, erratic nonphysical radar velocity values will appear in the solution. The problem requires some form of regularization based on the estimated errors in the data. Several approaches can be found in the recent literature for this purpose. Knoll and Clement (1999) adopted SVD, using only the 10 largest singular values. Buursink et al. (2002) used Tikhonov regularization after choosing the optimal number of layers on the basis of the Akaike information criterion, obtaining layers of 0.5-m thickness for their field site. An elegant and efficient solution for a similar problem was proposed by Lizarralde and Swift (1999), who extended the concept of Occam's inversion to vertical seismic profiles. Occam's inversion was developed originally for electromagnetic sounding data (Constable et al., 1987) and later adopted extensively for electrical resistivity tomography (LaBrecque et al., 1996). Occam's inversion for VRP finds the maximum value of the damping parameter {lambda} for which the objective function

[4]
is minimized to give

[5]
where

[6]
in which W is the diagonal matrix containing the reciprocals of data (travel time) error standard deviations and R is a roughness matrix, generally taken as a numerical approximation to the two-dimensional Laplacian operator. Constraint Eq. [5], which utilizes a {chi}2 statistic for the residuals in Eq. [6], is a direct consequence of the assumption that the travel time errors are normally distributed. In practical terms, Occam's inversion yields the smoothest (slowness) solution that is compatible with the error level in the data. Therefore, the error level must be estimated from field measurements, as discussed above. The model covariance matrix and the model resolution matrix can be computed as described by Alumbaugh and Newman (2000). The diagonal of the model covariance matrix contains the estimation variance for the radar velocity values. We utilized such variance to compute the confidence bands for radar velocity (and consequently moisture content; see Fig. 11) profiles, assuming a Gaussian probability distribution for the slowness, thus assuming that the mean value plus or minus two standard deviation values represent the 95% confidence interval.



View larger version (26K):
[in this window]
[in a new window]
  		Fig. 11. Comparison between vertical radar profiling (VRP)-derived moisture content profile for Borehole BB and the corresponding best result from Monte Carlo Richards' equation simulations: (a) 31 Mar. 2003; (b) 11 June 2003; (c) 26 Aug. 2003; (d) 3 Nov. 2003. The moisture content profile obtained with Occam's inversion is shown with a thick line, and relevant 95% confidence bands by thin lines. The best overall Monte Carlo simulation result is shown as a thick dashed line.

 
Petrophysical Models
A number of petrophysical models are potentially available to translate the bulk dielectric constant {kappa}r derived from Eq. [1] into moisture content estimates. The models can be broadly classified into two categories: (i) empirical relationships, such as the Topp et al. (1980) model that can be used without need for further soil characterization, and (ii) mixing models (Wharton et al., 1980; Roth et al., 1990; Knoll and Knight, 1994), for which site specific parameters need to be measured, such as solid matrix dielectric constant. For this work we adopted the Topp et al. (1980) relationship that is widely accepted in the TDR (and GPR) literature. We made a limited investigation on the use of the CRIM model, by assuming different solid matrix permittivity values. As discussed in the following sections, more important limitations in knowledge of the detailed subsurface structure makes the choice of the petrophysical relationship of lesser importance.

VRP Reflection Analysis
More information can be extracted from VRP data than only first arrivals. The energy arriving after the first break is often generated by reflections at lithological interfaces (Fig. 5) that are associated with sharp changes in the electromagnetic impedance and, ultimately, in permittivity values and radar velocities. It should be kept in mind that the tomographic inversion of the travel times via Occam's algorithm gives a smooth profile of the propagation velocity that cannot properly account for sudden changes in dielectric properties. However, the presence of sharp boundaries can be identified by the presence of reflections: the interfaces below the downhole antenna produce reflections of the direct signal (down-going wave-field), and the upward traveling wave is detected by the receiver antenna. This up-going wavefield can be easily identified in records since its slope is opposite to the down-going, direct wavefield (Fig. 5, 7) . Processing the records can enhance the reflections and produce an image of the interfaces in the subsoil, as in conventional vertical seismic profiling processing. The first step is to separate the up-going and down-going wavefields. This can be obtained with different approaches, such as f-k filtering (Yilmaz and Doherty, 2000). We based our processing on static corrections and median filter (Hardage, 2000). Raw data are corrected using a velocity model obtained from the picked travel times, which causes the down-going wavefield to be flattened. The median filter enhances the horizontal coherence, and the residual median filter extracts the up-going wavefield. A double opposite correction flattens the reflections, which then can be stacked using a corridor stacking. The stacked radargram gives an image of the interfaces in the subsurface; the identified sharp variations in electromagnetic impedance are interpreted as main lithological interfaces. In Fig. 8 we show the steps of reflection processing applied to the data collected 31 Mar. 2003, when the water table reached its lowest level at the Trecate site (10.68 m b.g.l.). This is the best condition for extracting reflection events from VRP data, since the deeper layers are also not fully saturated. Differences in lithological structure manifest themselves then as differences in the moisture content. Consequently, sharp contrasts in dielectric properties can exist, giving rise to reflections. Similar analyses conducted on data collected on other dates confirm the shallower reflection events (down to 6.5 m b.g.l.) but cannot show the deeper structure of the site. A similar reflection analysis was conducted on data from Borehole FW16 collected 14 Feb. 2003.



View larger version (76K):
[in this window]
[in a new window]
  		Fig. 8. Sequence of vertical radar profiling (VRP) processing to extract up-going reflections from data collected 31 Mar. 2003: (a) data after dewowing and static shift for instrument time-zero; (b) corrected for statics computed on first arrivals; (c) reflection enhancement via residual median filtering; and (d) reflections double corrected back with statics from first arrivals. From (d) the depth of the reflecting horizons can be read at the intersection between the corrected first arrival curve (black line) and each reflected event. Note that each reflection is represented by several peaks and troughs because of the wavelet of the radar signal. Uncertain reflectors are shown with a question mark.

 
In Fig. 8d, some reflected events can be clearly identified at 2, 6.5, and 8 m b.g.l.. The most notable reflections at 2 and 6.5 m are visible, even on the unprocessed data (Fig. 5). Note that each such event is characterized by more than one peak and trough, consistent with the characteristics of the radar wavelet (Fig. 3 and 4). Therefore, care must be taken (i) not to interpret multiple reflections when only one is present and conversely (ii) to identify reflections that are not well spaced in depth. For instance, a reflection at 9.5 m is likely to exist in Fig. 8, but the first break of the corresponding wavelet interferes with the energy train of the nearby reflection at 8.0 m. Another event is possible but difficult to confirm at 4 m. A fairly strong event at roughly 11 m b.g.l. seems to correlate with the position of the water table rather than to a stratigraphic boundary. The information derived from VRP reflections has been correlated with information from drilling logs (i.e., from the examination of rock cuttings during drilling itself). Such a process is known to be potentially affected by serious errors, especially since cuttings are brought to the surface by drilling fluid circulation, and its reliability depends heavily on the experience of the operators. In particular, while the existence of the main lithostratigraphic formations can be identified by drilling logs, the precise location of the stratigraphic boundaries can easily be misplaced. On the basis of this rationale, we corrected the drilling log results on the basis of the VRP reflection data (Fig. 9) . Note that in this process we discarded the possible reflection at the 4-m depth (Fig. 8) since this reflection has no supporting evidence from the drilling log information.



View larger version (23K):
[in this window]
[in a new window]
  		Fig. 9. Lithology at the Trecate site, derived from drilling logs and corrected via vertical radar profiling (VRP) reflection analysis for the two key boreholes considered in this study.

 
Infiltration Modeling
The result of the described VRP travel time inversion is a dataset of radar velocity profiles with depth at a number of points in time from December 2002 to November 2003. Such profiles are converted into moisture content profiles using Eq. [1] and the Topp et al. (1980) empirical model. This time–space dataset of moisture contents is an important piece of information that describes the dynamic hydrologic conditions of the deep vadose zone, in response to natural infiltration and water table oscillations. Hydrological models exist that can predict the moisture content as a function of space and time in response of the forcing boundary conditions. Such models are based on solutions of a nonlinear partial differential equation, known as the Richards' equation. Constitutive models that link the moisture content with water pressure or suction, and with the unsaturated hydraulic conductivity must be combined with the Richards' equation to apply the numerical solution to specific problems. Maybe the best known constitutive model is that of van Genuchten (1980), whose model depends on five parameters (two shape factors, {alpha}, n, a residual moisture content {theta}r, a saturated moisture content {theta}s, and the saturated hydraulic conductivity Ks). The identification of these five parameters for each formation in the subsurface defines uniquely the expected hydrological behavior of the system.

As discussed in our introduction, inversion of moisture content data derived from geophysical measurements to yield estimates of the governing hydrological parameters is becoming increasingly popular (Binley and Beven, 2003; Binley et al., 2002b, 2004; Kowalsky et al., 2004). We attempted a similar hydrogeophysical inversion using the VRP data at Trecate.

Hydrogeophysical inversion is similar to geophysical data inversion, in that it requires:

As a forward model, we adopted the widely available HYDRUS model (Simunek et al., 1998). Daily net rainfall estimates were used to define the surface boundary condition. Water table changes were simulated at the bottom of the soil column, honoring the 6-m yearly oscillations at the site. The transient simulations started on 1 Jan. 2002 after a 5-yr warm-up period used to bring the simulated soil column close to the site conditions. The initial condition on 1 Jan. 1997 was simulated using the steady-state simulator by Rockhold et al. (1997). The numerical grid consisted of 121 nodes at 10-cm spacing from the soil surface to 12 m b.g.l. Self-adjusting time steps were used, ranging from 0.0001 to 0.5 d, depending on convergence conditions of the model. Three lithologies were considered in the model as illustrated in Fig. 9. The adopted values for the van Genuchten parameters are shown in Table 1. The values for {alpha}, n, {theta}r, and {theta}s were derived from laboratory experiments on 5.08-cm (2 inch) Geoprobe cores extracted from depths ranging from 2 to 10 m b.g.l. The range in saturated hydraulic conductivity values was chosen to explore the expected intervals from the literature. An exhaustive Monte Carlo search strategy was used (i.e., parameters were chosen at random from a uniform probability distribution with selected parameter ranges). No correlation between parameter values was assumed to exist. We limited our Monte Carlo search to the saturated hydraulic conductivity values for the three materials. In this manner, we focused on the most important governing parameters that control the infiltration process, especially in relatively coarse materials. A complete Monte Carlo exploration of the full 3 x 5 parameter set would have required a considerably larger computational effort. A total of 20000 independent realizations were generated for each run. Less than 1% of the simulations did not converge because of the potentially high conductivity contrast across layers. For each Monte Carlo realization, the goodness of fit was computed with respect to VRP-derived moisture content profiles at all available time steps. The goodness of fit was quantified using an efficiency defined as 1 – {sigma}error2/{sigma}data2, where {sigma}error2 is the error of fit variance, and {sigma}data2 is the overall data variance.


View this table:
[in this window]
[in a new window]
  		Table 1. Ranges and values of van Genuchten parameters for the three material types.

 

   RESULTS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 MODELING
 RESULTS
CONCLUSIONS
 REFERENCES
 
The described procedure converted VRP profiles into quantitative moisture content profiles from which the hydraulic behavior of the unsaturated zone could be interpreted. Figure 10 shows the evolution of the moisture content according to VRP data for Borehole BB. Evident features are:



View larger version (33K):
[in this window]
[in a new window]
  		Fig. 10. Moisture content profiles at different times for Borehole BB, Trecate, derived from vertical radar profiling (VRP) via Occam's inversion and conversion from dielectric properties to moisture content.

 
The vadose zone at this site is highly dynamic, with considerable changes in the moisture profile at the same location between one survey and the next (about 2 wk). Some interesting results have been produced by comparing the Monte Carlo run results with the VRP-derived moisture content field data. Figure 11 shows a selection of simulated and measured moisture content profiles at the same time step. The simulated result is the one that best matches the moisture content at all time steps, and for both boreholes considered. Therefore, the match at each time step may not be the absolute best. In fact, Fig. 11 shows that it is possible to reproduce via Monte Carlo Richards' equation simulations the general pattern of the measured moisture content profiles, particularly their mean values over space at different time instants. However, the simulations differ from the measured profiles in the details. The simulated profiles show pronounced jumps in moisture content at the lithological boundaries, while the measured profiles are much smoother. Profiles at Borehole BB were generally reproduced more accurately than profiles at Borehole FW16 (not shown), for which the larger drilling size and gravel pack may have acted as a preferential pathway for faster infiltration. In fact, the average moisture content around FW16, as measured via VRP, was consistently lower than around Borehole BB. In Fig. 11 we show also the confidence bands relevant to the VRP-derived moisture content profiles. Note that such confidence bands are rather narrow, giving the impression that the accuracy of estimated moisture content profile is very good. It must be noted, though, that these confidence bands are those derived from the error estimates of the travel time, conditioned on the selected strategy for data inversion (i.e., the choice of the smoothest moisture content profile that honors the error level in the travel time data).

The choice of Occam's inversion to process VRP data may be the culprit for the discrepancies observed in Fig. 11 between the measured and simulated moisture content profiles. In fact, the choice of a smooth solution may be in contradiction with evidence that the moisture content profile is not smooth. The most convincing evidence for this is the existence of reflected events in the VRP data. Different VRP inversion approaches that do not rely on smoothness may prove to be more effective at reconstructing moisture content profiles at sites where such reflections are observed. The value of information in the VRP data, both in terms of (smoothed) moisture content profiles and in terms of reconstructed lithological boundaries, can be appreciated from Fig. 12 . Here we compare the observed and simulated moisture content profiles shown in Fig. 11a with a profile that can be derived from modeling in the absence of VRP data, by adopting the lithological boundaries defined by drilling logs and reasonable values of the vertical saturated hydraulic conductivity for material 3 (3.4 m d–1, as indicated by the analysis of pumping tests on site) and material 2 (0.5 m d–1).



View larger version (21K):
[in this window]
[in a new window]
  		Fig. 12. Comparison between the vertical radar profiling (VRP)-derived moisture content profile at Borehole BB (thick line and confidence band with thin lines), the corresponding best result from Monte Carlo Richards' equation simulations (thick dashed line), and a deterministic Richards' equation simulation based on the data available without VRP information (thick dot-dashed line). The deterministic simulation uses lithological boundaries derived from drilling logs, and the following hydraulic conductivity estimates: K1 = 10 cm d–1, K2 = 50 cm d–1, and K3 = 340 cm d–1 (value derived from in situ pumping test).

 
As discussed above, the goodness of fit between simulations and field data was quantified by an efficiency measure. The relevant dotty plots, representing overall efficiency of each simulation vs. the value of one governing parameter, are shown in Fig. 13 . The upper envelope of these dotty plots represents the best simulated results as a function of the parameter in abscissa. Flat upper envelopes mean that the parameter is not identified properly by the data and the model. This is the case for the saturated hydraulic conductivity of Material 1, the topsoil. Dotty plots that present pronounced peaks demonstrate that the corresponding parameter is well constrained. Figure 13 shows that the dotty plots for the hydraulic conductivity of Materials 2 and 3 have identifiable peaks. The same is true for the ratio of the hydraulic conductivity of the two materials. The values that yield the best efficiency (around 51%) are Ks = 0.61 m d–1 for Material 2 and Ks = 1.46 m d–1 for Material 3. Note that these values are fairly close to the only available estimate of the vertical hydraulic conductivity (under the water table) at the site (3.4 m d–1). These data are therefore in excellent agreement, considering that the saturated hydraulic conductivity often varies by orders of magnitude, even for the same geological formation.



View larger version (75K):
[in this window]
[in a new window]
  		Fig. 13. Dotty plots of the efficiency vs. (a) saturated hydraulic conductivity of Material 1 (topsoil); (b) saturated hydraulic conductivity of Material 2; (c) saturated hydraulic conductivity of Material 3; (d) ratio of saturated hydraulic conductivities of Materials 2 and 3. Note that while the hydraulic conductivity of the topsoil is poorly identified by the flat upper envelope of the dotty plot, the other three plots show distinct maxima.

 

   CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 MODELING
 RESULTS
 CONCLUSIONS
REFERENCES
 
Results from the Trecate site demonstrate the feasibility and value of the relatively simple and economic VRP approach for monitoring deep vadose zone dynamics. Data must be carefully processed, and special attention must be paid to obtaining good estimates of the travel time data error. This error estimate is essential for correct inversion of VRP travel times to yield radar velocity profiles. The vertical scale of measurement at this site, on the basis of the error analysis, is about 0.25 m. Also derived from the VRP data at Trecate was the location of lithological boundaries. This information can be extracted from processing and analysis of reflection events in the VRP radargrams. Both the moisture content profiles and the lithological boundaries derived from the VRP were used in a Monte Carlo hydrogeophysical inversion to yield estimates of the vertical saturated hydraulic conductivity of the key lithological units. From the comparison between simulations and field data, it appears that Occam's inversion may not be the best approach for reconstructing moisture content profiles at this site because the approach forces the results to be excessively smooth, especially in view of evidence that sharp contrasts exist that give rise to reflections.


   ACKNOWLEDGMENTS
 
The useful comments of two anonymous reviewers greatly improved the quality of the manuscript. We gratefully acknowledge Mark Rockhold, PNNL, USA, who granted us permission to use ss_infil. We also thank the Direzione Servizi Tecnici di Prevenzione, Regione Piemonte, for the meteorological data.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 MODELING
 RESULTS
 CONCLUSIONS
 REFERENCES
 




This article has been cited by other articles:


Home page
 Vadose Zone JHome page
R. Deiana, G. Cassiani, A. Villa, A. Bagliani, and V. Bruno
Calibration of a Vadose Zone Model Using Water Injection Monitored by GPR and Electrical Resistance Tomography
Vadose Zone J., February 25, 2008; 7(1): 215 - 226.
[Abstract] [Full Text] [PDF]

Home page
 Vadose Zone JHome page
M. C. Looms, K. H. Jensen, A. Binley, and L. Nielsen
Monitoring Unsaturated Flow and Transport Using Cross-Borehole Geophysical Methods
Vadose Zone J., February 25, 2008; 7(1): 227 - 237.
[Abstract] [Full Text] [PDF]

Home page
 Vadose Zone JHome page
M. C. Looms, A. Binley, K. H. Jensen, L. Nielsen, and T. M. Hansen
Identifying Unsaturated Hydraulic Parameters Using an Integrated Data Fusion Approach on Cross-Borehole Geophysical Data
Vadose Zone J., February 25, 2008; 7(1): 238 - 248.
[Abstract] [Full Text] [PDF]