Published in Vadose Zone Journal 3:1072-1081 (2004)
© 2004 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SPECIAL SECTION: HYDROGEOPHYSICS
Electromagnetic Inversion of GPR Signals and Subsequent Hydrodynamic Inversion to Estimate Effective Vadose Zone Hydraulic Properties
S. Lambota,*,
M. Antoinea,
I. van den Boschb,
E. C. Slobc and
M. Vancloostera
a Department of Environmental Sciences and Land Use Planning, Catholic University of Louvain, Croix du Sud 2, Box 2, B-1348 Louvain-la-Neuve, Belgium
b Microwave Laboratory, Catholic University of Louvain, Place du Levant 3, B-1348 Louvain-la-Neuve, Belgium
c Department of Geotechnology, Delft University of Technology, Mijnbouwstraat 120, 2628 RX Delft, The Netherlands
* Corresponding author (lambot{at}geru.ucl.ac.be)
Received 14 January 2004.
 |
ABSTRACT
|
|---|
We combine electromagnetic inversion of ground penetrating radar (GPR) signals with hydrodynamic inverse modeling to identify the effective soil hydraulic properties of a sand in laboratory conditions. Ground penetrating radar provides soil moisture time series that are subsequently used as input in the hydrodynamic inverse procedure. The technique relies on an ultrawide band (UWB) stepped frequency continuous wave (SFCW) radar combined with an off-ground monostatic transverse electromagnetic (TEM) horn antenna. Ground penetrating radar signal forward modeling is based on the exact solution of the three-dimensional Maxwell equations for describing free wave propagation and on linear systems in series and parallel for describing wave propagation in the antenna. Water flow in the sand is described by the one-dimensional Richards equation using the Mualem–van Genuchten parameterization. Both model inversions are formulated by the classical least-squares problem and are performed iteratively using advanced global optimization techniques. Compared with time domain reflectometry (TDR), results demonstrated the appropriateness of the GPR integrated approach to measure soil moisture remotely. In particular, the approach was found to be less sensitive to the inherent small-scale heterogeneities. Hydrodynamic inversion of soil moisture data led to hydraulic parameters agreeing reasonably well with direct measurements. The observed discrepancies were attributed to the different characterization scales and samples. The overall integrated approach offers great promise to map the effective hydraulic properties of the shallow subsurface at a high spatial resolution.
Abbreviations: GMCS, global multilevel coordinate search • GPR, ground penetrating radar • MVG, Mualem–van Genuchten model • NMS, Nelder–Mead simplex • SFCW, stepped frequency continuous wave • TDR, time domain reflectometry • TEM, transverse electromagnetic • UWB, ultrawide band • VNA, vector network analyzer
|
INTRODUCTION
|
|---|
CHARACTERIZATION OF vadose zone hydraulic properties remains a major challenge in the soil and hydrologic sciences. New appropriate methods are needed for many agricultural and environmental engineering applications involving water flow and solute transport modeling. It has become evident that modeling detailed spatial distributions of water and solutes in the heterogeneous subsurface requires extensive site characterization for determining the spatial distribution of the soil hydraulic properties, including the water retention curve and the hydraulic conductivity function (Yeh, 1998). Characterizing this variability with conventional laboratory or in situ downhole methods is invasive, and thus, time-consuming, costly, and subject to a high degree of uncertainty due to a lack of densely sampled in situ hydrological measurements. Given the limitations of these methods, intensive efforts have been undertaken to supplement the scarcity of hydrogeological data with densely sampled geophysical data (Beres and Haeni, 1991; Rubin et al., 1992; Hubbard et al., 1997; Hubbard and Rubin, 2000; Gloaguen et al., 2001).
In particular, GPR has been used increasingly to remotely identify soil stratigraphy (Davis and Annan, 1989; Kung and Lu, 1993; Boll et al., 1996), to locate the water table (Nakashima et al., 2001), to follow wetting front movement (Vellidis et al., 1990), to measure soil water content (Greaves et al., 1996; van Overmeeren et al., 1997; Weiler et al., 1998; Huisman et al., 2001; Rucker and Ferré, 2003; Grote et al., 2003), to assist in subsurface hydraulic parameter identification (Hubbard et al., 1997; Gloaguen et al., 2001), to assess soil salinity (al Hagrey and Müller, 2000), and to support the monitoring of contaminants (Brewster and Annan, 1994; Darayan et al., 1998; Yoder et al., 2001). An excellent review of GPR principles and history was given by Annan (2002), and reviews of its application for measuring soil water content in particular were given by Davis and Annan (2002) and Huisman et al. (2003).
Ground penetrating radar data can be utilized as information on the state of the flow system to infer the soil hydraulic properties using hydrodynamic inverse modeling methods (
imunek and van Genuchten, 1996; Durner et al., 1999; Romano and Santini, 1999; Hughson and Yeh, 2000; Hopmans et al., 2002). With these methods, state variables of the dynamic flow system are combined with a validated hydrodynamic model and an appropriate optimization algorithm to estimate the soil hydraulic parameters. The inversion is based on the optimization of an objective function that represents the error between the measured response of the soil system and the simulated response with the given model subject to a trial parameter set. Among existing hydrodynamic inverse modeling procedures, the method proposed by Lambot et al. (2002)(2004a) appears to be very promising to perform nondestructive hydraulic characterizations of the subsurface. Indeed, it requires only soil moisture time series gathered during a natural infiltration process as input. Soil moisture, in contrast to the conventionally used pressure head or flow data, is potentially obtainable from near-surface remote sensing using GPR.
The objective of this study is to combine the hydrodynamic inversion method of Lambot et al. (2002) with advanced GPR characterization techniques to identify the effective hydraulic properties of a sandy soil in laboratory conditions. The method provides the basis of a promising new integrated tool for laboratory and in situ characterization of the soil hydraulic properties. In addition, we compare GPR with direct conventional characterization methods and with TDR, which is already well established in the field of hydrodynamic inverse modeling. The main limitation of TDR is that it requires careful insertion of a probe into the soil. In this study, TDR is also used to monitor the soil electric conductivity. The electric conductivity, which is widely used for monitoring of contaminants, may be estimated from our GPR technique as well.
Following Lambot et al. (2004b)(2004c), the GPR system adopted in this study consists of an UWB–SFCW radar combined with a monostatic TEM horn antenna to be used off ground. This radar configuration is of practical interest because of its high mobility, but also because it enables realistic and efficient forward modeling of the radar–antenna–subsurface system (Lambot et al., 2003). As a result, the depth-dependent soil dielectric properties can be estimated from electromagnetic inversion of a single GPR signal (Lambot et al., 2004c). Traditional GPR techniques for measuring the depth-dependent soil water content usually require several measurements with two antennas in contact with the soil, or lowered into boreholes, to characterize a single profile (Nakashima et al., 2001; Garambois et al., 2002; Huisman et al., 2003). This severely restricts high-resolution and real-time mapping applications. To partly circumvent these limitations, Serbin and Or (2003) promoted the use of an off-ground monostatic horn antenna system and used it for identifying the water content of the surface soil skin (2 cm) from the surface reflection coefficient.
|
Ground Penetrating Radar System
|
|---|
Following Lambot et al. (2003)(2004b, 2004c), we used an UWB–SFCW radar combined with an off-ground monostatic traveling wave TEM horn antenna. A SFCW radar is a continuous wave system in which the carrier frequency is changed with a fixed step. A SFCW waveform is transmitted, sweeping a bandwidth at discrete frequency values. A TEM horn antenna consists of a pair of triangular conductors forming a V structure guiding essentially a TEM mode between its antenna plates (i.e., the electric and magnetic fields are transverse to the propagation direction).
This radar configuration has several advantages over commonly used time domain technologies and ground-based bistatic antenna systems. The major value of a SFCW radar stems from the potential for controlling an UWB, which is needed for better depth resolution and to obtain more information from the subsurface. Commercially available time domain GPR systems generally have a bandwidth limited to <1 GHz. Moreover, the signal/noise ratio and the dynamic range of SFCW radars is better than for time domain systems since the mean radiated power is much higher and each sampled frequency is an independent measurement.
The high directivity (focused beaming) of the TEM horn yields high spatial resolution data and enables off-ground measurements, which ensures a high degree of mobility. Also, the ground can be situated in the Fraunhofer region (far-field) of the antenna, which can then be modeled simply and accurately as a point source and receiver, instead of a spatially distributed source and receiver (Hansen and Bailin, 1959; Balanis, 1997). The monostatic mode permits efficient and realistic forward modeling of the radar–antenna–subsurface system, since in this mode the received signal propagates mostly along the antenna axial direction. As a result, antenna modeling does not require knowledge of the source radiation pattern. Also, the stochastic horizontal variability of the dielectric properties inherently encountered in environmental systems is expected to play a negligible role, such that the ground can be modeled accurately by means of a horizontally multilayered configuration, for which closed form analytical expressions can be derived by solution of the system of Maxwell's equations. Finally, in monostatic mode, there is no coupling (direct air wave) between the antennas. In bistatic mode of operation, direct signals are the main source of noise.
|
Inversion of the GPR Signal
|
|---|
Antenna Equation in the Frequency Domain
In this study, a monostatic UWB–SFCW radar system was set up using a vector network analyzer (VNA). As a result, the GPR signal to be modeled consists of the frequency-dependent complex ratio S11(
) between the returned signal and the emitted signal, being the angular frequency. The VNA reference plane where S11() is actually measured was established at the connection between the antenna feed point and the VNA cable.
Lambot et al. (2004c) proposed the block diagram depicted in Fig. 1
as the antenna model. The model consists in a simple linear system composed of elementary model components in series and parallel, all characterized by their own frequency response function accounting for specific electromagnetic phenomena. The resulting transfer function relating S11() measured by the VNA to the frequency response Gxx
() of the multilayered medium (air–subsurface system) is expressed in the frequency domain by
 | [1] |
where Y() and X() are, respectively, the received and emitted signals at the VNA reference plane; Hi(), Ht(), Hr(), and Hf() are the complex return loss, transmitting, receiving, and feedback loss transfer functions of the antenna, respectively; and Gxx
is the transfer function of the air–subsurface system modeled as a multilayered medium. Equation [1] is referred to as the antenna equation in the frequency domain. In the time domain, the scalar products in Eq. [1] correspond to time convolution integrals.
View larger version (17K):
[in this window]
[in a new window]
|
Fig. 1. Block diagram representing the vector network analyzer–antenna–multilayered medium system modeled as linear systems in series and parallel.
|
|
Due to variations in the impedance between the antenna feed point, antenna aperture, and air, multiple wave reflections (ringing) occur within the antenna. The return loss transfer function Hi() represents the part of this ringing, measured at the reference plane, that is independent of the backscattered electromagnetic field Gxx. The transmitting and receiving transfer functions describe the antenna gain and phase delay between the measurement point and the source and receiver virtual point. The positive feedback loop with transfer function Hf(), similarly to Hi(), accounts for impedance variations between the antenna feed point, antenna aperture, and air, which cause a part of the backscattered field to be reflected again toward the subsurface. This leads to multiple wave reflections between the antenna and the subsurface. We refer the reader to Lambot et al. (2004c) for additional details on estimating these characteristic antenna response functions.
Green's Function
The air–subsurface system is modeled as a three-dimensional stratified dielectric and conductive planar system consisting of N horizontal layers (vertical in our study) separated by N – 1 interfaces parallel to the 
plane, as illustrated in Fig. 2
. The nth layer is homogeneous and characterized by a dielectric permittivity 
n, electric conductivity 
n, and thickness hn. The magnetic permeability µ is considered constant and equal to the free space permeability µ0. The emitting part of the TEM horn antenna is approximated by an infinitesimal horizontal x-directed electric dipole [i.e., the second subscript x in Gxx, whereas the receiving part of the antenna is emulated by recording the horizontal x-directed component of the backscattered electric field [first subscript x in Gxx]. The point source is located at the origin S of the coordinate system. The use of a three-dimensional model is essential to account for the spherical divergence (geometric spreading) in wave propagation.
The theoretical basis for GPR wave propagation is given by Maxwell's equations. Green's function (i.e., the solution of Maxwell's equations) for electromagnetic waves propagating in multilayered media is well known (Arbel and Felsen, 1963; Tai, 1994; Michalski and Mosig, 1997; Peterson et al., 1998). Following the approach of Slob and Fokkema (2002) and Lambot et al. (2004c), the analytic expression for Green's function in the spectral domain (two-dimensional spatial Fourier domain) for our specific source and receiver configuration is found to be
 | [2] |
where the subscript n equals 1 and denotes here the first interface and first layer; RnTM and RnTE are the transverse magnetic and transverse electric global reflection coefficients accounting for all reflections and multiples from subsurface interfaces, respectively; 
n is the vertical wave-number defined as n = 
where kp is a spectral domain transform parameter, 
n = jµn, 
n = n + jn, and j = –1.
The global TM-mode and TE-mode reflection coefficients at interface n (n = 1, ..., N – 1) are given by
 | [3] |
 | [4] |
 | [5] |
 | [6] |
where rnTM and rnTE denote the local-plane wave TM and TE mode reflection coefficients, respectively, at interface n. These expressions are in a recursive form. The recursion is initiated by the observation that there are no upgoing waves from the lower (left) half-space, such that RN–1TM = rN–1TM and RN–1TE = Rn–1TE. The local reflection coefficients determine that part of the electromagnetic wave reflected at a dielectric interface, while the other part is transmitted.
The transformation of Eq. [2] from the spectral domain to the spatial domain is performed by employing the two-dimensional Fourier inverse transformation:
 | [7] |
which reduces to a single integral in view of the invariance of the dielectric properties along the x and y coordinates and the monostatic mode configuration. The integral can be computed with standard integration routines using Gaussian-like quadrature. Green's function (Eq. [7]) represents the x component of the backscattered electric field 
at the virtual source point of the antenna for a unit strength x-directed electric source.
Since in this study we are only interested in measuring the soil water content, which is directly related to the dielectric permittivity, we consider an inverse modeling analysis in which the frequency-dependent dielectric parameters are replaced by frequency-independent effective parameters. Lambot et al. (2004c) previously showed that the dielectric permittivity of sand, in contrast to the electric conductivity, is almost frequency independent in the frequency range considered. This simplification reduces the number of parameters to be estimated in the electromagnetic inversion procedure.
|
Objective Function
|
|---|
Subsurface parameter identification by inverse modeling is a nonlinear optimization problem that consists of finding the parameter vector b = [n,n,hn]T, n = 1···N, such that an objective function 
(b) is minimized. For the particular case where no prior information about the parameters is used and assuming that observation errors are normally distributed with mean zero and covariance matrix C, independent, and homoskedastic, the maximum likelihood theory reduces to a weighted least squares problem. The objective function to be minimized is defined accordingly as follows (Bard, 1974; Press, 1992):
 | [8] |
where Gxx = Gxx
and Gxx = Gxx
are vectors containing, respectively, the observed and simulated response functions. Since the response function is a complex function, the difference between observed and modeled data is expressed by the amplitude of the errors in the complex plane.
Objective function Eq. [8] indirectly relates the response function of the multilayered medium to its constitutive parameters. However, as for most electromagnetic inverse problems, this function is highly nonlinear and characterized by oscillatory behavior and a multitude of local minima. This complex topography necessitates the use of a robust global optimization algorithm. Following the approach of Lambot et al. (2002)(2004b), we used the global multilevel coordinate search (GMCS) algorithm (Huyer and Neumaier, 1999) combined sequentially with the classical Nelder–Mead simplex algorithm (NMS) (Lagarias et al., 1998).
|
Hydrodynamic Modeling and Inversion
|
|---|
The governing flow model relates the subsurface hydraulic properties to the response of the system (i.e., the soil moisture time series), which constitutes simultaneously in this study the output of the GPR signal electromagnetic inversion and the input of the hydrodynamic inversion. One-dimensional vertical transient water flow in homogeneous isotropic rigid porous media is described by the Richards equation, expressed here in terms of pressure head (Jury et al., 1996):
 | [9] |
where h is the time- and depth-dependent pressure head, C(h) = 
(h)/h is the differential water capacity with (h) being the water retention curve and the volumetric water content, K() is the hydraulic conductivity, t is time, and z is depth taken positive downward. The relations (h) and K() are the constitutive functions that describe the soil hydraulic properties. In this study, Eq. [9] is solved subject to the following initial and boundary conditions:
 | [10a] |
 | [10b] |
 | [10c] |
that prescribe at the upper boundary a zero flux condition, and at the lower boundary a variable positive pressure head, with L being the soil profile length.
The nonhysteretic unimodal Mualem–van Genuchten model (MVG) (Mualem, 1976; van Genuchten, 1980) is used in this study to describe the soil hydraulic properties. The water retention curve is given by
 | [11] |
where r and s are the residual and saturated water contents, respectively; 
and n are curve shape parameters inversely related to the air-entry value and the width of the pore size distribution, respectively; and m is restricted by the Mualem condition m = 1 – 1/n, with n > 1. The hydraulic conductivity relationship is given by
 | [12] |
where Ks is the saturated hydraulic conductivity, Se = ( – r)/(s – r) is effective saturation, and 
is a factor that accounts for pore connectivity or pore tortuosity. A high value of corresponds to a low pore connectivity or a high pore tortuosity.
Differential Eq. [9] is solved numerically by finite differences using a MATLAB (The MathWorks, 1999) implementation of the WAVE model (Vanclooster et al., 1996). In this code, a fully implicit scheme is used to advance the solution in time with automatically adjusting time steps on the basis of a mass conservative convergence criterion. Instantaneous relative errors in the mass balance were <10–4. We discretized the spatial flow domain into 93 equidistant elements each representing 2 cm of the 186-cm-deep sand profile.
The total number of unknown hydraulic parameters in the MVG model is six, that is, r, s, , and n for the water retention curve plus Ks and for the unsaturated hydraulic conductivity. s has a clear physical significance and can easily be obtained by direct measurements. r is usually defined as the residual water content corresponding to a value of h 
–
. Generally, this parameter is regarded as an empirical parameter and can be fixed either to a value that yields the best fit to the experimental water retention data (Kool et al., 1985), or to zero (Nimmo, 1991; Fuentes et al., 1992). Following Lambot et al. (2002), we measured s directly and fixed r to zero, thus reducing the unknown parameter vector to b = [,n,Ks,]T.
Assuming measurement errors to be Gaussian distributed with mean zero and covariance matrix C, and assuming errors to be independent and homoskedastic, the hydrodynamic inversion is formulated by a weighted least-squares problem, for which the objective function is defined by (Bard, 1974; Press, 1992)
 | [13] |
where * = *(z,t) and = (z,t,b) are the vectors (size n x 1) containing the observed and simulated water content time series, respectively; e is the vector of residuals; and n is the total number of observations. Covariance matrix C is a diagonal matrix whose elements are identical and equal to the variance of the measurement errors. Objective function Eq. [13] is minimized using the same optimization procedure as for the electromagnetic inversion process, that is, using GMCS–NMS (Lambot et al., 2002).
|
MATERIALS AND METHODS
|
|---|
A schematic representation of the laboratory experimental setup is depicted in Fig. 3
. A tank (1.90 m high, 1 m wide, and 0.1 m thick) was filled with sand. To avoid air entrapment, the sand was packed under water saturated conditions. This resulted in a slight layering similar to natural deposits. The side facing the antenna was made of a 1-cm-thick transparent PVC sheet, with known dielectric properties. To control the boundary conditions of the electromagnetic model, the opposite side of the tank consisted of a metal sheet, which was assumed to reflect perfectly electromagnetic waves. Consequently, extraneous laboratory materials behind this metal sheet did not influence the measured backscattered GPR signal. The antenna was situated about 40 cm from the tank, with measurements performed with the antenna situated at two different depths, namely, z1 and z2. The measurement scale of the GPR system depends on the average operating frequency (f = 2 GHz in this study) and the distance between the antenna and the sand (h = 0.40 m). The measurement scale is expressed by the Fresnel zone whose radius Fr can be computed as
 | [14] |
where c = 3 x 108 m s–1 is the free space propagation speed of electromagnetic waves. In the hydrodynamic model, we considered a linearly weighted averaging scheme for the GPR measure of water content, with the weighting coefficient equal to 1 at the center of the Fresnel zone and to 0 at its limits.
We implemented an UWB–SFCW radar system using a VNA (ZVRE, Rohde & Schwarz, Munich, Germany) with an excellent dynamic range (>130 dB). The antenna system consisted of a linear polarized double ridged broadband TEM horn (BBHA 9120 A, Schwarzbeck Mess-Elektronik, Schönau, Germany). The antenna was 22 cm long with a 14 x 24 cm2 aperture area, and the antenna weight was 1.3 kg. Its nominal frequency range was 0.8 to 5 GHz. The antenna was connected to the reflection port of the VNA via an N type 50-
impedance coaxial cable. We calibrated the VNA at the connection between the antenna feed point and the cable using a 50- OSM (Open, Short, Match) series of a high precision standard calibration kit (ZVZ21-N, Rohde & Schwarz). This established a reference calibration plane to which the frequency-dependent complex ratio S11 between the returned signal and the emitted signal was measured. Parameter S11 was measured sequentially at 301 stepped operating frequencies for the range of 0.8 to 3.2 GHz, with a frequency step of 8 MHz.
Additionally, the tank was equipped with two TDR probes inserted horizontally into the sand layer at the two monitoring depths, at about 35 cm from the center of the tank to avoid interference with GPR. Time domain reflectometry measurements were performed with a Tektronix 1502C metallic TDR cable tester (Tektronix Inc., Beaverton, OR). The probes were constructed in our laboratory and consisted of three stainless-steel rods (9.3 cm long and 0.5-cm diameter), with a 2.5-cm rod spacing. The probes were connected to the Tektronix via 2-m-long coaxial cables and BNC connectors. The Tektronix was connected to a computer with an RS232 interface to numerically acquire the TDR signal. The soil dielectric permittivity and electric conductivity were acquired using the WinTDR 98 program (Soil Physics Group, Utah State University, Logan, UT; Or et al., 1998). In the hydrodynamic model, the TDR measurement scale was assumed to correspond to the model discretization (i.e., 2 cm).
Starting from complete saturation, the sand was allowed to drain by progressively decreasing the bottom pressure head to zero. Soil moisture time series were monitored using both TDR and GPR. The water content was derived from the dielectric permittivity using the sand specific relations presented in Fig. 4
. These relations were determined by Lambot et al. (2004c) for laboratory conditions by the same GPR and TDR techniques as used here. The slight difference between the TDR and GPR calibration curves may be due to the frequency dependence of the dielectric permittivity. The higher frequencies of GPR lead to higher relaxation phenomena, which is expressed by a lower apparent dielectric permittivity for the same water content, particularly at the higher water content levels. The standard deviation of the error in the water content prediction using GPR was 0.0070 m3 m–3, which reflects the high accuracy of the method. For TDR, the observed standard deviation was 0.0032 m3 m–3.
Direct measurements of the sand hydraulic properties were made with a sandbox apparatus on 100-cm3 samples to provide the drying water retention curve in the suction range 0 to 80 cm, and a constant head apparatus to give the saturated conductivity Ks on a 1000-cm3 sample (e.g., see Klute and Dirksen, 1986). The parameter in the MVG model was set to the commonly used value of 0.5. The hydraulic parameters of the sand are reported in Table 1 (Direct).
View this table:
[in this window]
[in a new window]
|
Table 1. Mualem–van Genuchten parameters of the sand for the synthetic (inversion of error-free and error-contaminated data) and real experiments (inversion of GPR and TDR data).
|
|
|
RESULTS AND DISCUSSION
|
|---|
Hydrodynamic Inversion of Synthetic Data
Using synthetic experiments with error-free and error-contaminated data we first demonstrate that measured soil moisture time series at the two specified depths for our specific laboratory conditions leads to a unique and stable solution for the hydrodynamic inverse problem. The boundary conditions used for the synthetic experiments were the same as for the real experiment. The soil properties (true parameters) used to generate the synthetic data were those of the actual sand determined by conventional direct laboratory methods. Inversions were performed in a large parameter space given by [0.010 < < 0.050 cm–1, 1.10 < n < 15, 1 x 10–2 < Ks < 1 x 101 cm min–1, 0 < < 5].
Results are presented in Table 1. We can observe that the parameters obtained by inversion of the error-free data correspond exactly to the true parameters. This demonstrates that enough information is contained, at least theoretically, within the soil moisture time series to ensure a unique solution. Synthetic water content data were subsequently corrupted with normally distributed errors with a mean of 0.01 m3 m–3 and a standard deviation of 70 x 10–4 m3 m–3 (the standard deviation for the GPR data shown in Fig. 4). Results show that the impact of these errors on the inversely estimated parameters is acceptable, which indicates the good stability properties of the considered inverse problem. Further details on the well-posedness of the adopted inversion scheme are reported in Lambot et al. (2002).
GPR Signal Inversion
Lambot et al. (2004b) demonstrated that the electromagnetic inverse problem dealt with in this study is well-posed. Figure 5 represents the measured and modeled Green's functions taken at an arbitrary time during the drainage experiment. Results are similar for the other measurements. For convenience, data are presented in both the frequency and time domains. The time domain shows two main reflections. The first, arising between 3 and 4 ns, stems from the PVC sheet and the sand surface, while the second, arising between 5 and 6 ns, is due to the metal sheet. The water content of the soil directly affects the electromagnetic wave propagation time between these two reflectors. The oscillating behavior of the time domain signal stems from the limited frequency range considered for the inverse Fourier transform (i.e., 0.8–3.2 GHz). It is worth noting that estimation of the dielectric permittivity from only the two-way travel time observed in the time domain would not be as accurate as the inverse modeling procedure. Indeed, it can be demonstrated numerically from the exact solution of Maxwell's equations that the positions of the peaks or the zero-crossing points in the time domain does not correspond exactly to the positions of the interfaces. Moreover, it is not possible to derive the electric conductivity from the travel time analysis.
View larger version (20K):
[in this window]
[in a new window]
|
Fig. 5. Measured and modeled Green's functions represented in the (a) frequency domain and (b) time domain.
|
|
The measured and modeled frequency domain Green's functions agree relatively well. Given the experimental setup in relation to the electromagnetic model assumptions, we attribute the origin of the observed discrepancies mainly to (i) the frequency dependence of the dielectric properties of the sand, which were not accounted for in the present study; (ii) the ambiguous noise due to extraneous scattering by objects present in the laboratory; (iii) the finite size of the tank (the approximation of an infinite medium may be incorrect, particularly at the lower frequencies); and (iv) errors in determining the antenna transfer functions. However, we note that the effect of these different sources of error on the water content should be negligible since, as illustrated by the time domain representation of the Green's functions, the propagation time is very well modeled. The propagation time is directly related to the phase of the frequency domain complex Green's function, which is closely reproduced, except for frequencies where the amplitude of the Green's function is low or tends to go to zero. In this case, small errors on the real and imaginary parts of the Green's function can result in large errors in the phase. This is an artifact of representing the complex data in terms of modulus and phase.
Soil moisture time series derived from GPR signal inversion are compared in Fig. 6
with time series obtained from TDR measurements. Results for GPR are very consistent regarding the drainage experiment. The water content at depth z1 first decreases, followed by the water content at depth z2, slowly reaching a state of hydrostatic equilibrium with a water table at depth zb (Fig. 3). For TDR, the results are less consistent. After the first 50 min of drainage, the water content pertaining to depth z2 reached lower values than the water content at depth z1. This is attributed mainly to variability in the sand properties at the TDR measurement scale. Our packing of the tank led to a slight layering in the sand, which may have influenced local soil moisture and water flow dynamics. Because of its larger measurement scale, GPR is less susceptible to this inherent local variability. Additionally, systematic TDR measurement errors are not to be excluded. Indeed, different travel time analysis methods can lead to differences as high as 0.04 m3 m–3 for the same waveform (Huisman et al., 2002). Also, travel time analysis parameters are often fixed by users who then visually test whether or not the analysis is performed correctly. However, it is difficult to obtain a set of travel time analysis parameters that do not need to be changed during processing of TDR measurements for a wide range of soil water contents. In our study we did not change those parameters.
Water Content Time Series Inversion
Modeled soil moisture time series obtained by hydrodynamic inversion are presented in Fig. 6. The GPR time series were well reproduced by the hydrodynamic model. In contrast, the TDR data were not well reproduced because of important modeling and/or measurement errors, as explained in the above section. The hydraulic parameters obtained by model inversion are depicted in Table 1 for both TDR and GPR and compared with the parameters obtained using conventional direct laboratory methods. A difference of one order of magnitude occurs between the directly determined saturated conductivity and the GPR derived estimate. This type of error in the saturated conductivity is commonly observed for direct measurements (Schaap and Leij, 2000). Bruckler et al. (2002) also observed such discrepancies, and even larger, in estimated hydraulic conductivities using inversion of infiltration data. The prediction from TDR data was significantly overestimated (two orders of magnitude). This is due to the combined effect of a large modeling error in the prediction of the water content time series, and high correlation between Ks and (r = 0.9816). The effect of a large Ks on the unsaturated hydraulic conductivity is compensated by a large , which denotes substantial tortuosity in soil structure. We attribute the differences in the hydraulic properties determined by the three methods mainly to the different samples perhaps having different structural properties (bulk density).
The corresponding water retention and hydraulic conductivity functions are presented in Fig. 7
. The inversely estimated curves are all very characteristic of a sandy soil, but with some important differences. In particular, the saturated water contents for GPR and TDR are higher than for the direct measurements. Most of the discrepancies stem from measurement and calibration errors inherent to the different characterization techniques, but also from the different measurement scales. For anisotropic and heterogeneous porous media as inherently encountered in the environment, effective soil physical properties depend on the measurement scale (Yeh et al., 1985; Khaleel et al., 2002). This finding reflects the variability in the sand bulk density among the different samples.
View larger version (15K):
[in this window]
[in a new window]
|
Fig. 7. (a) Water retention curve and (b) hydraulic conductivity function of the sand determined directly, by inversion of GPR data, and by inversion of TDR data.
|
|
Electric Conductivity
Figure 8
represents relations between the electric conductivity and dielectric permittivity for the GPR and TDR measurements pertaining to depth z1. Notice that, as for TDR, the GPR derived electric conductivities are very consistent with the different dielectric values or water contents. This means that the GPR technique used is appropriate for measuring this variable as well, and hence offers great promise for monitoring solute transport. The relations observed for TDR and GPR are slightly different. This can be explained by the different operating frequencies of GPR and TDR, the different measurement scales, and the different measurement locations involved. These relations depend very much on soil structure and texture (e.g., see Rhoades et al., 1976). Lambot et al. (2004c) obtained far less satisfactory results since measurements were performed on different mixtures of the sand.
|
SUMMARY AND CONCLUSIONS
|
|---|
In this study, we combined the hydrodynamic inverse modeling procedure of Lambot et al. (2002) with the advanced GPR integrated approach of Lambot et al. (2004c) to identify the effective soil hydraulic properties of a sandy soil under laboratory conditions. Electromagnetic inversion of the GPR signal provides soil moisture data, which are subsequently inverted with a hydrodynamic model to estimate the hydraulic parameters of interest.
Ground penetrating radar and TDR measurements were performed along a tank filled with a sandy soil subject to a natural drainage from initially saturated conditions. An advantage of GPR compared with TDR is that it operates at a larger scale, and consequently is less influenced by inherent small-scale heterogeneities. Ground penetrating radar hence leads to soil moisture time series that are more consistent with the homogeneity assumption in the hydrodynamic model. Reducing modeling errors in hydrodynamic inversion procedures is important since they usually constitute the main source of error in hydraulic parameter identification.
The inversely estimated hydraulic properties agreed reasonably well with properties obtained using traditional direct methods. The observed discrepancies were attributed mostly to the different characterization scales and samples involved. The estimated saturated conductivity was one order of magnitude higher for GPR. The GPR derived water retention curve was slightly shifted compared with the directly determined curve, but still was very typical of a sand.
We believe that the proposed integrated approach constitutes a valuable alternative to the commonly used methods and represents the basis of a new integrated tool with great promise for laboratory and in situ characterizations of the soil hydraulic properties. The approach has been found applicable to the identification of the effective properties of large-scale laboratory soil samples. Moreover, as for TDR, the GPR approach permits measurements of the electric conductivity, and therefore also offers promise for solute transport modeling. Our intention is to focus next on inversion of GPR data for field conditions, which is complicated due to unknown depths of the interfaces and topographic roughness of the air–soil interface.
|
ACKNOWLEDGMENTS
|
|---|
This work was supported by a research grant from the "Fonds pour la Formation à la Recherche dans l'Industrie et dans l'Agriculture" (Belgium), the "Fonds National de la Recherche Scientifique" (Belgium), the Delft University of Technology (The Netherlands), the Royal Military Academy (Belgium), and the Catholic University of Louvain (Belgium).
|
REFERENCES
|
|---|
- al Hagrey, S., and C. Müller. 2000. GPR study of pore water content and salinity in sand. Geophys. Prospect. 48:63–85.
- Annan, A. 2002. GPR—History, trends, and future developments. Subsurf. Sens. Technol. Appl. 3:253–270.
- Arbel, E., and L.B. Felsen. 1963. Theory of radiation from sources in anisotropic media. Part 1: General sources in stratified media. In E. Jordan (ed.) Electromagnetic theory and antennas. Part 1. Macmillan, New York.
- Balanis, C.A. 1997. Antenna theory: Analysis and design. John Wiley and Sons, New York.
- Bard, Y. 1974. Nonlinear parameter estimation. Academic Press, New York.
- Beres, M., and F. Haeni. 1991. Application of ground penetrating radar methods in hydrogeologic studies. Ground Water 29:375–386.[ISI]
- Boll, J., R. van Rijn, K. Weiler, T. Steenhuis, J. Daliparthy, and S. Herbert. 1996. Using ground penetrating radar to detect layers in a sandy field soil. Geoderma 70:117–132.[ISI]
- Brewster, M., and A. Annan. 1994. Ground penetrating radar monitoring of a controlled DNAPL release: 200 MHz radar. Geophysics 59:1211–1221.[ISI]
- Bruckler, L., P. Bertuzzi, R. Angulo-Jaramillo, and S. Ruy. 2002. Testing an infiltration method for estimating soil hydraulic properties in the laboratory. Soil Sci. Soc. Am. J. 66:384–395.[Abstract/Free Full Text]
- Darayan, S., C. Liu, L. Shen, and D. Shattuck. 1998. Measurement of electrical properties of contaminated soils. Geophys. Prospect. 46:477–488.
- Davis, J., and A. Annan. 1989. Ground penetrating radar for high resolution mapping of soil and rock stratigraphy. Geophys. Prospect. 37:531–551.
- Davis, J.L., and A.P. Annan. 2002. Ground penetrating radar to measure soil water content. p. 446–463. In J.H. Dane and G.C. Topp. Methods of soil analysis. Part 4. SSSA Book Ser. 5. SSSA, Madison, WI.
- Durner, W., B. Schultze, and T. Zurmül. 1999. State-of-the-art in inverse modeling of inflow/outflow experiments. p. 661–681. In M.Th. van Genuchten et al. (ed.) Proceedings of the International Workshop on Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media. University of California Press, Riverside, CA.
- Fuentes, C., R. Haverkamp, and J.-Y. Parlange. 1992. Parameter constraints on closed-form soil–water relationships. J. Hydrol. (Amsterdam) 134:117–142.
- Garambois, S., P. Sénéchal, and H. Perroud. 2002. On the use of combined geophysical methods to assess water content and water conductivity of near-surface formations. J. Hydrol. (Amsterdam) 259:32–48.
- Gloaguen, E., M. Couteau, D. Marcotte, and R. Chapuis. 2001. Estimation of hydraulic conductivity of an unconfined aquifer using cokriging of GPR and hydrostratigraphic data. J. Appl. Geophys. 47:135–152.
- Greaves, R., D. Lesmes, J. Lee, and M. Toksov. 1996. Velocity variations and water content estimated from multi-offset, ground-penetrating radar. Geophysics 61:683–695.[ISI]
- Grote, K., S. Hubbard, and Y. Rubin. 2003. Field-scale estimation of volumetric water content using GPR ground wave techniques. Water Resour. Res. 39(11):1321. doi:10.1029/2003WR002045
- Hansen, R., and L. Bailin. 1959. A new method of near field analysis. IRE Trans. Antennas Propag. AP-7:S458–S467.
- Hopmans, J.W., J. imunek, N. Romano, and W. Durner. 2002. Inverse methods. p. 963–1008. In J.H. Dane and G.C. Topp (ed.) Methods of soil analysis. Part 4. SSSA Book Ser. 5. SSSA, Madison, WI.
- Hubbard, S., and Y. Rubin. 2000. Hydrogeological parameter estimation using geological data: A review of selected techniques. J. Contam. Hydrol. 45:3–34.
- Hubbard, S., Y. Rubin, and E. Majer. 1997. Ground-penetrating-radar-assisted saturation and permeability estimation in bimodal systems. Water Resour. Res. 33:971–990.
- Hughson, D., and T.-C. Yeh. 2000. An inverse model for three-dimensional flow in variably saturated porous media. Water Resour. Res. 36:829–839.
- Huisman, J., S. Hubbard, J. Redman, and A. Annan. 2003. Measuring soil water content with ground penetrating radar: A review. Available at www.vadosezonejournal.org. Vadose Zone J. 2:476–491.[Abstract/Free Full Text]
- Huisman, J., C. Sperl, W. Bouten, and J. Verstraten. 2001. Soil water content measurements at different scales: Accuracy of time domain reflectometry and ground penetrating radar. J. Hydrol. (Amsterdam) 245:48–58.
- Huisman, J., A. Weerts, T. Heimovaara, and W. Bouten. 2002. Comparison of travel time analysis and inverse modeling for soil water content determination with time domain reflectometry. Water Resour. Res. 38:1224. doil:10.1029/2001WR000259.
- Huyer, W., and A. Neumaier. 1999. Global optimization by multilevel coordinate search. J. Global Optim. 14:331–355.
- Jury, W., W. Gardner, and W. Gardner. 1996. Soil physics. 5th ed. John Wiley and Sons, New York.
- Khaleel, R., T.-C. Yeh, and Z. Lu. 2002. Upscaled flow and transport properties for heterogeneous unsaturated media. Water Resour. Resear. 38(5):1053. doi:10.1029/2000WR000072
- Klute, A., and C. Dirksen. 1986. Conductivities and diffusivities of unsaturated soils. p 687–734. In A. Klute (ed.) Methods of soil analysis. Part 1. 2nd ed. Agron. Monogr. 9. ASA and SSSA, Madison, WI.
- Kool, J., J. Parker, and M.Th. van Genuchten. 1985. Determining soil hydraulic properties from one-step outflow experiments by parameter estimation: Theory and numerical studies. Soil Sci. Soc. Am. J. 49:1348–1354.[Abstract/Free Full Text]
- Kung, K.-J., and Z.-B. Lu. 1993. Using ground penetrating radar to detect layers of discontinuous dielectric constant. Soil Sci. Soc. Am. J. 57:335–340.[Abstract/Free Full Text]
- Lagarias, J., J. Reeds, M. Wright, and P. Wright. 1998. Convergence properties of the Nelder-Mead Simplex method in low dimensions. Siam J. Optim. 9:112–147.
- Lambot, S., F. Hupet, M. Javaux, and M. Vanclooster. 2004a. Laboratory evaluation of a hydrodynamic inverse modeling method based on water content data. Water Resour. Res. 40(3):W03506. doi:10.1029/2003WR002641
- Lambot, S., M. Javaux, F. Hupet, and M. Vanclooster. 2002. A global multilevel coordinate search procedure for estimating the unsaturated soil hydraulic properties. Water Resour. Res. 38(11):1224. doi:10.1029/2001WR001224
- Lambot, S., E.C. Slob, I. van den Bosch, B. Stockbroeckx, B. Scheers, and M. Vanclooster. 2003. GPR design and modeling for identifying the shallow subsurface dielectric properties. p. 130–135. In A. Yarovoy (ed.) Proceedings of the 2nd International Workshop on Advanced Ground Penetrating Radar. Delft University of Technology, Delft, The Netherlands.
- Lambot, S., E. Slob, I. van den Bosch, B. Stockbroeckx, B. Scheers, and M. Vanclooster. 2004b. Estimating soil electric properties from monostatic ground-penetrating radar signal inversion in the frequency domain. Water Resour. Res. 40:W04205. doi:10.1029/2003WR002095.
- Lambot, S., E. Slob, I. van den Bosch, B. Stockbroeckx, and M. Vanclooster. 2004c. Modeling of ground penetrating radar for accurate characterization of subsurface dielectric properties. IEEE Trans. Geosci. Remote Sens. TGRS-00347-2003. (In press.)
- MathWorks. 1999. MATLAB. Version 5.3. MathWorks, Natick, MA.
- Michalski, K., and J. Mosig. 1997. Multilayered media Green's functions in integral equation formulations. IEEE Trans. Antennas Propag. 45:508–519.
- Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522.
- Nakashima, Y., H. Zhou, and M. Sato. 2001. Estimation of groundwater level by GPR in an area with multiple ambiguous reflections. J. Appl. Geophys. 47:241–249.
- Nimmo, J. 1991. Comment on the treatment of residual water content in "A consistent set of parametric models for the two-phase flow of immiscible fluids in the subsurface" by L. Luckner et al. Water Resour. Res. 27:661–662.
- Or, D., B. Fisher, R. Hubsher, and J. Wraith. 1998. WinTDR98 v4.0—Users guide. Plant, Soils, & Biometeorology, Utah State Univ., Logan, UT.
- Peterson, A.F., S.L. Ray, and R. Mittra. 1998. Computational methods for electromagnetics. IEEE–Oxford University Press, New York.
- Press, W. 1992. Numerical recipes in FORTRAN: The art of scientific computing. 2nd ed. Cambridge University Press, New York.
- Rhoades, J., P. Raats, and R. Prather. 1976. Effects of liquid-phase electrical conductivity, water content, and surface conductivity on bulk soil electrical conductivity. Soil Sci. Soc. Am. J. 40:651–655.[Abstract/Free Full Text]
- Romano, N., and A. Santini. 1999. Determining soil hydraulic functions from evaporation experiments by a parameter estimation approach: Experimental verifications and numerical studies. Water Resour. Res. 35:3343–3359.
- Rubin, Y., G. Mavko, and J. Harris. 1992. Mapping permeability in heterogeneous aquifers using hydrologic and seismic data. Water Resour. Res. 28:1809–1816.
- Rucker, D., and T. Ferré. 2003. Near-surface water content estimation with borehole ground penetrating radar using critically refracted waves. Available at www.vadosezonejournal.org. Vadose Zone J. 2:247–252.[Abstract/Free Full Text]
- Schaap, M., and F. Leij. 2000. Improved prediction of unsaturated hydraulic conductivity with the Mualem–van Genuchten model. Soil Sci. Soc. Am. J. 64:843–851.[Abstract/Free Full Text]
- Serbin, G., and D. Or. 2003. Near-surface water content measurements using horn antenna radar: Methodology and overview. Available at www.vadosezonejournal.org. Vadose Zone J. 2:500–510.[Abstract/Free Full Text]
- imunek, J., and M. van Genuchten. 1996. Estimating unsaturated soil hydraulic properties from tension disc infiltrometer data by numerical inversion. Water Resour. Res. 32:2683–2696.
- Slob, E., and J. Fokkema. 2002. Coupling effects of two electric dipoles on an interface. Radio Sci. 37:1073. doi:10.1029/2001RS2529
- Tai, C.-T. 1994. Dyadic Green function in electromagnetic theory. IEEE Press, New York.
- Vanclooster, M., P. Viaene, J. Diels, and K. Christiaens. 1996. WAVE, a mathematical model for simulating water and agrochemicals in the soil and vadose environment. Release 2.1. Institute for Land and Water Management, K.U.Leuven, Belgium.
- van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.[Abstract/Free Full Text]
- van Overmeeren, R., S. Sariowan, and J. Gehrels. 1997. Ground penetrating radar for determining volumetric soil water content: Results of comparative measurements at two test sites. J. Hydrol. (Amsterdam) 197:316–338.
- Vellidis, G., M. Smith, D. Thomas, and L. Asmussen. 1990. Detecting wetting front movement in a sandy soil with ground penetrating radar. Trans. ASAE 33:1867–1874.
TOP
ABSTRACT
INTRODUCTION
