Published in Vadose Zone Journal 3:1128-1145 (2004)
© 2004 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SPECIAL SECTION: HYDROGEOPHYSICS
Obtaining the Spatial Distribution of Water Content along a TDR Probe Using the SCEM-UA Bayesian Inverse Modeling Scheme
Timo J. Heimovaaraa,*,
Johan A. Huismanb,c,
Jasper A. Vrugtb and
Willem Boutenb
a Royal Haskoning, P.O. Box 8520, 3009 AM Rotterdam, The Netherlands
b Institute for Biodiversity and Ecosystem Dynamics (IBED), University of Amsterdam, Nieuwe Achtergracht 166, 1018 WV Amsterdam, The Netherlands
c Institute for Landscape Ecology and Resource Management, Justus Liebig University, Heinrich-Buff-Ring 26-32, D-35392 Giessen, Germany
* Corresponding author (T.Heimovaara{at}RoyalHaskoning.com)
Received 13 May 2004.
 |
ABSTRACT
|
|---|
Time domain reflectometry (TDR) has become one of the standard methods for the measurement of the temporal and spatial distribution of water saturation in soils. Current waveform analysis methodology gives a measurement of the average water content along the length of the TDR probe. Close inspection of TDR waveforms shows that heterogeneity in water content along the probe can be seen in the TDR waveform. We present a comprehensive approach to TDR waveform analysis that gives a quantitative estimate of the dielectric permittivity profile along the length of the probe and, therefore, the distribution of water content. The approach is based on the combination of a multisection scatter function model for the TDR measurement system with the shuffled complex evolution Metropolis algorithm (SCEM-UA). This combined approach allows for the estimation of the 40 parameters in the transmission line model using a series of simple calibration measurements. The proof of concept is given with measurements in a layered system consisting of air and water. Finally, TDR waveforms from layered soil samples were analyzed to estimate the distribution of the water content along the length of the probe. Results show that the proposed method provides much more reproducible results than obtained with the traditional travel time method. Because the proposed method can be fully automated, it increases the applicability of the TDR method, especially in applications where detailed (real-time) data are required on heterogeneous infiltration.
Abbreviations: MSSF, multi-section scatter function • SCEM-UA, shuffled complex evolution Metropolis algorithm • SR, scale-reduction [score] • TDR, time domain reflectometry
|
INTRODUCTION
|
|---|
KNOWLEDGE OF THE spatial and temporal variability of water saturation in soils is important to obtain improved estimates of water flow (and its dissolved components) through the vadose zone. Because of its high accuracy and potential for automated measurement, TDR has become one of the standard methods to measure the spatial and temporal variability of water contents in laboratory soil cores and experimental field plots. The method is based on the measurement of a reflected voltage-wave from a probe installed in the soil. Detailed analysis of this reflected waveform enables the estimation of the bulk-electrical conductivity and dielectric permittivity of the soil, which in porous media primarily depends on the water content. Traditionally, waveforms are analyzed such that an average water content estimate of the soil along the TDR probe is obtained. In this work we show that measured TDR waveforms contain sufficient information to estimate water content profiles along TDR probes, corresponding with a scale of cubic centimeters (10–6 m3) to cubic decimeters (10–3 m3).
Two different approaches to obtaining soil water content estimates from measured TDR waveforms have found widespread use. The first estimates the dielectric permittivity of the soil (and thus the water content) by detecting characteristic reflections in the waveform (Topp et al., 1980; Heimovaara and Bouten, 1990; Hook et al., 1992; Robinson et al., 2003). In particular, the travel time is used to estimate the dielectric permittivity, while the amplitude changes in the waveform are used to assess the bulk electrical conductivity. This method is relatively simple to implement and can be used for all known TDR probes (Heimovaara and Bouten, 1990; Baker and Allmaras, 1990; Jones et al., 2002; Robinson et al., 2003). However, in the presence of layered soils and porous media with a relatively high electrical conductivity, the method exhibits difficulties with the correct interpretation of the waveforms, often leading to erroneous water content estimates (i.e., Schaap et al., 2003). Coating the wave-guides with an insulating material can solve some of the problems related to high electrical conductivity (Ferré et al., 1996). Nevertheless, the intrinsic uncertainty of this method requires a substantial commitment of human resources to find and correct outliers in the results due to erroneous (automated) analyses. This is a problem for fully automated measurement and analysis. The second, perhaps more sophisticated method is based on the modeling of the waveform using the physical properties of the TDR system (Giese and Tiemann, 1975; Heimovaara, 1994; Heimovaara et al., 1996; Friel and Or, 1999; Feng et al., 1999; Schlaeger et al., 2001; Lin, 2003a, 2003b). In this approach, the dielectric permittivity and bulk electrical conductivity are model parameters estimated by fitting against measured waveforms. Although the method requires significant computer resources, the number of outliers due to erroneous analyses is substantially reduced (Huisman et al., 2002).
We follow the latter approach, using deterministic modeling of the TDR system in the frequency domain. This inverse modeling method originated in the 1970s (Giese and Tiemann, 1975; Clarkson et al., 1977) and has been improved in the 1990s (Heimovaara, 1994; Heimovaara et al., 1996). Recently, major steps forward have been made so that it is possible to explicitly account for the multiple sections in the measurement setup with the so-called multi-section scatter function (MSSF) (Feng et al., 1999; Lin, 2003a, 2003b).
While this inverse modeling approach results in a more objective analysis of TDR waveforms, with less commitment of human resources, the successful implementation of this method is critically dependent on the availability of optimization algorithms, which can reliably solve for high-dimensional parameter estimation problems. The successful application of computerized optimization techniques becomes more complicated as the number of parameters increases; this is because of the presence of local minima in the parameter space and the possibility that parameters are highly correlated (Duan et al., 1992). For example, in its simplest form, the MSSF model already contains eight unknown parameters. An increased number of sections to more accurately describe the underlying TDR–soil system leads to an even larger number of model parameters that needs to be estimated by fitting against measured TDR waveforms. Fortunately, significant advances have recently been made in the field of global nonlinear optimization. The shuffled complex evolution Metropolis global optimization algorithm (SCEM-UA) is a general-purpose code that reliably locates the global optimum in the parameter space for a high-dimensional optimization problem. In addition, the probability density functions of all parameters in the vicinity of the global optimum are also obtained (Vrugt et al., 2003).
Our objective is to demonstrate that the combined use of the MSSF model with the SCEM-UA algorithm facilitates the estimation of water content profiles along a standard TDR probe. With this general framework, it is possible to improve estimates of the spatial variability of water saturation on the scale of a TDR measurement. To achieve this goal we used the SCEM-UA algorithm for both the calibration of the measurement system as well as for the inverse modeling of the dielectric properties of the soil. The insight in the probability density of the optimized parameters provides the opportunity to assess the uncertainty of our experimental findings. The applicability of the proposed method is quite general because all measurements were performed with standard TDR instrumentation suitable for application in the field.
We present a comprehensive approach for an advanced form of TDR measurement interpretation. After summarizing the theoretical background we describe our approach, consisting of a calibration procedure in three steps and subsequently the measurement procedure. The measurement procedure is demonstrated with two cases, measurements in a layered system with air and water and measurements on layered soil samples.
|
THEORY
|
|---|
Modeling TDR waveforms in the frequency domain is based on the assumption that the measured time-domain waveform, r(t), is a convolution of the incident or input signal, v0(t), and a system function, s(t), describing how the system under test influences the input signal. The convolution theorem states that a convolution in the time domain is a multiplication in the frequency domain:
 | [1] |
in which R(f) signifies the Fourier transform of the (measured) time-domain waveform, r(t); V0(f) is the Fourier transform of the incident signal generated by the cable tester, v0(t); and S(f) denotes the Fourier transform of the system function, s(t). Knowledge of V0(f) and S(f) allows R(f) to be calculated, from which the corresponding time domain waveform is obtained (Heimovaara, 1994; Friel and Or, 1999).
To calculate r(t) or R(f), models are required to describe S(f), the so-called scatter function, V0(f), the input-function. Also, because S(f) depends on the complex frequency-dependent dielectric relaxation properties of the materials within the transmission line, a model is also required for the frequency-dependent dielectric relaxation. The derivation of the MSSF model for S(f) is presented in brief in the Appendix, as it has been thoroughly discussed previously (Heimovaara, 1994; Friel and Or, 1999; Feng et al., 1999).
The Input Signal, V0(f)
Heimovaara (2001) used an analytical function to model the input signal in the time domain:
 | [2] |
in which erf is the error function, t denotes time, 
is a parameter signifying the inverse of the rise time, and t0 is the position where the input signal starts rise. This model gives a symmetrical step function for which the rise time and position of the step can easily be modified. To apply this input function in the frequency domain, the function is transformed with a fast Fourier transform using the same approach as was used by Heimovaara (1994) and Friel and Or (1999).
The use of Eq. [2] to model the input signal has several advantages. First, we no longer require a measured input signal with identical time resolution as the TDR waveform measurements (Heimovaara, 1994; Friel and Or, 1999; Feng et al., 1999; Weerts et al., 2001). Instead, we can use the MSSF to model the effect of the measurement system on the incident signal before it enters the probe. Second, we do not have to measure the complete waveform as we can simulate the complete waveform and compare the simulated waveform to the window of measurement. As a result the resolution of TDR measurements can be adapted to the complexity of the measured signal to capture all the required detail.
The Frequency-Dependent Complex Dielectric Permittivity
Several models can be used to describe the frequency-dependent complex dielectric permittivity of a material (Heimovaara et al., 1994; Lin, 2003b). In this study we assumed that all materials in our transmission line exhibit a single relaxation in the frequency range of our measurement. We describe this dielectric relaxation with the Debeye relaxation function (Hasted, 1973):
 | [3] |
where 
s is the relative static permittivity, H is the relative high-frequency permittivity (value at the highest frequency in the measurement range), and frel signifies the relaxation frequency, defined as the frequency at which the permittivity equals (s + H)/2.
Waveform Calculation Procedure
To simulate TDR waveforms with the MSSF model we require a time-axis spanning the complete TDR time domain (which is from just before the incident signal is picked up by the sampling device to at least the end of the measured domain). The time axis is obtained in two steps from the time axis of a measured waveform. The first step is to add points to the front of the time axis so the first point falls before –2 ns. The second step is to add points to the end of the time axis so the number of points becomes a power of 2 for an efficient use of fast Fourier transforms (Press et al., 1986). The waveform calculation procedure is outlined in Heimovaara (1994) and Friel and Or (1999). To prevent aliasing in the Fourier and inverse Fourier transform it is particularly important that the time step in the time axis used for the calculation of the waveforms is <75 ps. This corresponds to a sampling frequency that is >6.6 GHz.
SCEM-UA
The SCEM-UA was recently developed by Vrugt et al. (2003). This algorithm is an adaptive evolutionary Monte Carlo Markov Chain method inspired by the SCE-UA global optimization algorithm of Duan et al. (1992) and combines the strengths of the Metropolis algorithm (Metropolis et al., 1953), controlled random search (Price, 1987), competitive evolution (Holland, 1975), and complex shuffling (Duan et al., 1992) to obtain an efficient estimate of the most optimal parameter set, and its underlying posterior distribution, within a single optimization run. The SCEM-UA algorithm is based on a Bayesian inference scheme (Thiemann et al., 2001) in which the a priori probability density distribution of the model parameters is updated to a posterior distribution. The posterior distribution is given with p(
|Y), in which is the vector of model parameters and Y denotes a vector of measurements. The prior information usually consists of fixing lower and upper boundary values for each of the parameters in , thereby creating the feasible parameter space, and specifying a uniform density over this hypercube. A uniform prior density distribution with Guassian error residuals leads to the following form of the posterior density:
 | [4] |
where Yj is the jth term of m measurements and 
j denotes the corresponding model prediction using the parameter vector . Notice, that this form of the posterior density minimizes the sum of squared errors between measurements and model predictions.
To generate samples from Eq. [4], the SCEM-UA algorithm starts by randomly selecting an initial population of points from the uniform prior parameter distributions using Latin hypercube sampling. For each of these points the posterior density is calculated using Eq. [4]. The population is then partitioned into a number of complexes using the calculated posterior densities. From each of the complexes a Markov chain is evolved using the Metropolis–Hastings algorithm (Metropolis et al., 1953; Hastings, 1970). After a preset number of iterations, the information contained in each of the sequences is shared by shuffling all complexes and restarting the procedure. These algorithmic steps are performed until all sequences have converged to a limiting posterior distribution. The advantage of the SCEM-UA approach is that the method conducts an efficient search of the feasible parameter space for high-dimensional problems.
To test for convergence of the parallel sequences to a limiting posterior distribution, we used two convergence diagnostics. First the scale-reduction (SR) score developed by Gelman and Rubin (1992) is required to be <1.2, as was proposed by Vrugt et al. (2003). However, practical experience with the calibration of the MSSF model lead us to suggest that this measure is not sufficient to test for convergence alone. For example, in some calibration runs significant better parameter combinations were continued to be found in the high probability density region of the parameter space, even though the SR value already indicated convergence to the posterior target distribution. Therefore, we also required that either the maximum difference in the last 25 generated SR scores was <0.001 or that the SR score was <1.03. Additionally, we also set a maximum to the number of model evaluations (iterations), ranging between 25000 and 100000, depending on the complexity of the optimization problem. A large number of iterations are required for high-dimensional parameter estimation problems, and poorly known prior distributions.
Estimation of Water Content
Topp et al. (1980) presented a calibration equation relating the apparent dielectric permittivity Ka to the volumetric soil water content (w):
 | [5] |
The apparent dielectric permittivity is closely related to the real part of the dielectric permittivity in the upper range of the TDR measurement bandwidth around 1GHz (Heimovaara et al., 1996, 2002; Lin, 2003b). We assume that Ka 
Re[(1GHz)].
|
MATERIALS AND METHODS
|
|---|
General Approach
To illustrate the power and applicability of the combined MSSF model and SCEM-UA inverse modeling framework we performed a number of measurements, which are summarized in Table 1. The first measurements in Step 1 to 3 are used to calibrate the TDR measurement system. After calibration, we performed an experiment in which we took measurements of a TDR probe partly submerged in water to different depths. The last experiments are on layered-soil samples and were included to illustrate the usefulness and applicability of the proposed inverse modeling method for estimating the water-saturation distribution along a TDR probe. Each set of measurements was analyzed with the SCEM-UA algorithm to obtain the posterior probability distribution of the parameters. For completeness, for each of the measurements, Table 1 also lists the number of parameters that were estimated with the SCEM-UA algorithm. In the next paragraphs, we present a condensed description of our measurement setup, summarize the details of how we compared measured waveforms with calculated waveforms, describe the calibration approach, and finally discuss the different experiments that were performed.
TDR Setup and Probes
Measurements were performed using a Tektronix (Beaverton, OR) 1502B cable tester connected to a PC through the Tektronix SP232 serial interface. We used a triple-wire TDR probe as described in Heimovaara (1993). The length of the probe was approximately 0.1 m, and the inner wire was about 0.007 m shorter than the two outer wires. The triple-wire probe was connected to the cable tester with RG 58 C/U-type coaxial cable with a length of 3 m.
Time domain reflectometry waveforms were obtained with the automated TDR measurement system developed by Heimovaara and Bouten (1990). This system stores waveform data together with all cable-tester settings. Two types of waveforms were stored, standard waveforms with 251 points that are normally used for the traditional time-domain analysis and waveforms with 1024 points that were used for the frequency-domain analysis described by Heimovaara (1994) and Heimovaara et al. (1996).
Calculation of Difference between Modeled and Measured Waveforms
The density criterion specified in Eq. [4] is calculated using the difference (residuals) between measured and modeled waveforms. Figure 1
illustrates waveforms with 251 and 1024 points from measurement of a layered sample. The enlargement shows the full range of the 251 points. It is obvious that the waveform with 251 points captures more detail than the waveform with 1024 points. The measurement software automatically optimizes the cable-tester settings to capture the first reflection with the highest resolution across 251 points. The lower resolution of the 1024-point waveform results in a smoothing of sharp reflections (marked by the arrows in Fig. 1). During inverse modeling, we used all measured values from the 251-data point waveform augmented with points from the 1024-waveform that lie beyond the 251-waveform time range. In this way we saved a considerable amount of measurement time, as this approach does not require a measurement of the whole waveform with the highest possible resolution.
View larger version (15K):
[in this window]
[in a new window]
|
Fig. 1. High-resolution 251-point waveform (circles) and low-resolution 1024-point waveform (crosses). The arrows in the inset mark waveform features missed by the 1024-point waveform.
|
|
Calibration Measurements (Step 1 to 3)
The first calibration step (Step 1 in Table 1) is needed to determine the two parameters ( and t0) for the analytical input function in Eq.[2]. We measured a waveform with a 50-
load (supplied with the Tektronix 1502 series cable testers) while making sure that the initial step in the waveform was sampled. In the second step we measured waveforms for three open-ended RG 58 C/U type coaxial cables with different lengths (0.99, 1.34, and 3.22 m, respectively). These measurements were used to obtain the transmission line parameters of the cable tester, the coaxial connectors and the coaxial cables (Eq. [A3] to [A9]). In the third calibration step we performed measurements in air and water with a known temperature. These measurements were used to obtain the transmission line parameters of the triple wire TDR probe (Eq. [A11] and [A12]).
Inverse Modeling for Calibration
Figure 2
presents a schematic overview of the TDR measurement system, which is used throughout the remainder of this paper. The system with the triple-wire probe is modeled with six transmission-line sections. The system with the open-ended cable lengths is modeled with four sections (the two-probe sections are not required). Figure 2 also indicates the numbering convention used in this paper (similar to the one used in Feng et al., 1999). The rationale for counting from the probe back toward the cable tester lies in the algorithm presented in the Appendix. Sections 3 to 6 (the cable, connector and internal part of the cable tester) are modeled using Eq. [A3] to [A9], as these are fully coaxial sections. The two-probe sections are modeled with Eq. [A11] to [A13], as the probe is not fully coaxial but consists of three wires.
View larger version (11K):
[in this window]
[in a new window]
|
Fig. 2. Overview of the multisection transmission line connected to the cable tester. The numbers indicate the transmission-line section numbering convention adopted for this paper following the approach of Feng et al. (1999). Section 4 is the transmission-line section inside the cable tester, Section 3 is the coaxial cable connecting the probe to the CT, Section 2 is the probe to cable interface, and finally Section 1 is the probe.
|
|
For the coaxial transmission line sections, Eq. [A10] shows that Z0 is calculated from the ratio of b and a. Knowledge of two of these parameters allows the third to be calculated. To reduce the number of parameters in the MSSF model, we chose to set the value for b to 7.25 mm to allow Z0 to be determined by inverse modeling. The value of 7.25 mm corresponds to the dimensions of a RG213 high-quality coaxial cable. We note that this value is chosen completely arbitrarily, but in all circumstances will be compensated by the value of a. In the absence of compelling information on the internal length of the transmission-line sections in the cable tester, we set the value of s of the internal cable-tester sections and the connector to 5. Preliminary calibration runs with the SCEM-UA algorithm, demonstrated increased convergence rates to a limiting posterior distribution when fixing this parameter because this parameter is highly correlated with the length of the internal sections.
Table 2 gives a summary of the parameters that must be estimated by calibration. In total we require 40 parameters before we are able to use the presented approach to measure the frequency-dependent dielectric properties of a (soil) sample. The parameters were identified in three steps (Table 1).
In the first step the two parameters for the input function were obtained by using SCEM-UA to fit the input function to the measurement of the 50- load. The mean values from the posterior probability distributions are used in the next steps.
In the second step the transmission-line parameters of the internal cable-tester sections, the connector, and the three cables were estimated. The parameters of transmission-line sections three to six in Eq. [A3] to [A9] were estimated by fitting the three open-ended cable measurements simultaneously with a four-section (fully coaxial) model using the parameters for the analytical input function from the first step. In this step we also assumed that s, of the coaxial cable was equal to 5 and that the parameters for the three coaxial cables are the same except for Ll and Z0. As a result we required the simultaneous estimation of 32 parameters (7 x 3 = 21 for Sections 4, 5, and 6, and 5 + 3 x 2 = 11 for Section 3).
The third step was to obtain the parameters for the measurement setup with the triple-wire probe by fitting the model to measurements of the probe in air and in water. To have the best possible calibration of our measurement setup we determined the distribution of all 38 transmission-line parameters (shown in Table 2) by using the posterior distributions from the second calibration step to estimate the prior parameter distribution for Sections 3 to 6. We also decided to optimize the electrical conductivity of our water sample, even though we have directly measured values. We did this so we could use the optimized electrical conductivity to verify the model of the transmission line. In addition, not requiring a separate measurement of the electrical conductivity of the water makes the calibration one step simpler. Table 3 lists the prior and posterior ranges of the parameters in the MSSF model after calibration of the TDR measurement system; the prior range is marked with the subscript "ini". The allowable search range lies between "Min" and "Max".
Experiment 1: Measurement of the Air–Water Boundary along the Probe
To assess the accuracy of measuring layers of different permittivity with a single probe by inverse modeling, we measured a series of waveforms of a probe inserted into water to different depths. To accurately determine the insertion-depth of the probe we mounted the probe on a high-precision workbench, which allowed us to measure the depth of insertion with an accuracy of approximately 1 µm. A series of waveforms was recorded, starting with the probe fully in air, then inserted into water to different depths until fully submerged, and finally with the probe fully in air again.
The waveforms of this layered air–water system were analyzed using inverse modeling with the SCEM-UA algorithm. We modeled this system by assuming the probe section of our transmission-line system to consist of two sections with a combined length equal to the length of the probe. In this particular case, there were seven transmission-line sections. The first section is that part of the probe that is inserted in water, while the second section is the part that is surrounded with air. The only parameter that was optimized was the relative length of the air section.
Experiment 2: Measurement of the Permittivity Distribution in Packed Layered Samples
To verify whether the recorded TDR waveforms contain sufficient information to estimate the spatial distribution of water saturation along a TDR probe, we created a packed sample of wet loam with a gravimetric water content of 0.26 g g–1 and air-dry fine sand with a gravimetric water content of 0.004 g g–1. Samples were carefully packed around the wires of the probe fixed in stainless-steel rings of 0.05-m diameter and 0.10-m height. Two samples were packed, both samples with approximately 0.05 m loam and 0.05 m sand. The first sample (loam–sand) was packed with the layer of loam closest to the probe head followed by a layer of sand at the end of the probe. The second sample was packed with a layering sequence of sand and loam. TDR measurements were performed on both samples and then the gravimetric water content of the samples was measured by oven drying. The waveforms obtained during Exp. 2 were analyzed using inverse modeling with the SCEM-UA algorithm. We used two approaches. In the first approach we assumed that the sample consists of two different layers. In the second approach we assumed that the sample consists of four or eight layers.
In the first approach, the relative lengths of the layers, and the associated Debeye parameters and 
DC for each layer were optimized by fitting against the corresponding measured waveform. Because of the small jitter in measured waveforms we also optimized the length of the cable and the characteristic impedance of the cable to have the most accurate description of the waveform. For the two-layer model we required one parameter for the relative length. The result of this is an 11-parameter optimization problem for the two-layer model; all other parameters were fixed at the mean values of the posterior distributions obtained from the calibration steps.
The prior parameter ranges for the two-layer model were obtained from an automated qualitative interpretation of the waveform, in which we used the amplitude of the waveform to estimate the maximum and minimum dielectric permittivity with Eq. [A11]. During optimization, the parameters were allowed to vary in a wider range between 1 and 300 for s and H, between –10 and 1 for log(DC), between 0 and 10.5 for log(frel), and between 0 and 1 for the length fraction of the different probe sections. We used the high value for s because the static permittivity in soils can reach very high values due to the so-called Maxwell–Wagner relaxation of bound-water layers or double layers (Hilhorst, 1998).
In the second waveform analysis approach we increased the number of layers in our model to 4 and 8 sequentially. We used the results from the two- and four-layer optimizations to determine the prior parameter ranges for the four- and eight-layer optimizations. The prior ranges were estimated using the mean value with a dispersion of three standard deviations. The cable length and characteristic impedance for the cable were fixed to the mean value of the two-layer optimization. The number of parameters optimized was 16 for the four-layer model and 32 for the eight-layer model.
|
RESULTS AND DISCUSSION
|
|---|
Calibration Step 1: The Input Waveform
The first step of the calibration procedure was to determine the parameters of the input function (Eq. [2]) to be used to model the input signal of the cable tester. Figure 3
shows the calibrated fit of the input signal, both in the frequency domain and in the time domain. The dots (heavy line in the frequency charts) show the measured signal; the lines are obtained from the input model. The corresponding parameter values derived with SCEM-UA are given in Table 4.
View larger version (20K):
[in this window]
[in a new window]
|
Fig. 3. Measured and simulated input signals based on optimized parameters. The top chart shows a small section of the input signal in the time domain. The middle and bottom charts show the magnitude and phase of the complex input signal as a function of the complete frequency range. Solid lines show the results for the model; the open circles are points in the measured waveform.
|
|
View this table:
[in this window]
[in a new window]
|
Table 4. Calibrated parameters for the input function of the Tektronix 1502B cable tester used for the TDR measurements.
|
|
These parameters indicate that the travel time between the voltage-sampling device and the zero position of the cable tester is approximately 1.6206 ns. The rise time of the input signal is 1/ and is therefore 156 ps. This value agrees very well with the specifications of the cable tester (Tektronix, 1988). The differences between the measured and simulated input signal can be attributed to sections in the transmission line within the cable tester that have an impedance value different from the assumed 50- characteristic impedance. For example, a small mismatch is seen at 0 ns due to the connector between the 50- load and the cable tester.
Calibration Step 2: Cable Tester and Cable Parameters Using Open-End Cables
Figure 4
shows the measured and simulated waveforms of the three open-end cables after all 32 parameters are optimized. In general, the multisection transmission-line model can accurately describe the measurements. The addition, Eq. [A9] to the MSSF model makes it possible to describe the dispersion due to the skin effect, marked in Fig. 4. We also see that the model fails to describe the interface between the cable tester and the cable at travel times close to 0 s. Nonetheless, this is not a very disturbing problem, as the window of dielectric measurements focuses on the primary and secondary reflections. These are described accurately so we consider the model to be quite accurate, if not correct. The fit of the longest cable is not perfect. This could be caused by the fact that we assumed the transmission-line parameters for all three cables to be identical except for the cable length (Ll) and characteristic impedance (Z0).
View larger version (17K):
[in this window]
[in a new window]
|
Fig. 4. Waveforms from open-end cables used to optimize the parameters for the internal transmission-line section of the cable tester and the RG58-C/U coaxial cable. Open circles indicate measurements; solid lines are model fits.
|
|
The SCEM-UA derived mean values and standard deviations of the parameters for the three-cable sections are given in Table 5. Most of the parameters in the MSSF model are well determined by calibration against measured TDR waveforms. Hence, the posterior parameter ranges are narrow for most of the model parameters, indicating that the TDR waveform contains sufficient information to identify 32 model parameters. Notice that the optimized value for Z0 differs from one cable to another. Seemingly, to obtain an optimal fit to the measurements, parameters need to be optimized for each measurement setup individually.
View this table:
[in this window]
[in a new window]
|
Table 5. Mean values and standard deviations (italics) from the optimal probability distributions obtained with the SCEM-UA algorithm for the three open-end cables. The parameters L and Z0 were optimized for each cable separately, and s was set to five for all transmission-line sections. All other parameters were assumed to be the same for all three cables.
|
|
Calibration Step 3: Complete TDR Measurement System
Figure 5
shows the graphical comparison of modeled and measured waveforms of the probe in air and in water. It is clear that the proposed inverse modeling approach using the MSSF model and the SCEM-UA algorithm allows for an accurate description of the measured waveforms.
View larger version (14K):
[in this window]
[in a new window]
|
Fig. 5. Measured and modeled waveforms in air and water after calibration. The solid lines are the results from the four-section model; open circles are the sampled points of the waveform sections. The top chart is of the high-resolution 251-point measurements with optimized cable-tester resolution. The bottom chart shows the 1024-point measurements with fixed cable-tester resolution.
|
|
Table 3 summarizes the mean values and the standard deviations of all 38 parameters obtained with the SCEM-UA algorithm for the water and air measurements. The small values for the standard deviation indicate a small uncertainty in the optimized parameter values. The electrical conductivity of the water was measured with a portable field-sensor to be 0.0563 S m–1. The mean value for the optimized electrical conductivity was 0.0569 S m–1 with a 95% confidence interval between 0.0566 and 0.0571 S m–1. The relative difference between the two methods is <1% and falls well within the accuracy of the field sensor. The implication of this result is that the model used to simulate our TDR measurement setup is physically correct, which increases the faith in the correctness of the estimated values of the other parameters. Another important implication from a practical viewpoint is that no additional electrical conductivity measurement is required for the calibration of the measurement setup.
Although the number of parameters that was determined by inverse modeling is very large (38), the results clearly indicate that the measured TDR waveforms contain enough information to precisely determine the transmission-line parameters. Moreover, the parameters obtained make sense from a physical point of view, especially considering the relatively large prior ranges of the model parameters in which the SCEM-UA algorithm was allowed to search.
Even though the SCEM-UA algorithm is a global search algorithm, good prior estimates are essential for a fast convergence to a limiting distribution. Indeed, although not further demonstrated here, using information from the three open-end cables to define the prior parameter ranges significantly decreased the number of model evaluations (about 25000) needed for the SCEM-UA algorithm to reach convergence.
Experiment 1: Layered System of Air and Water
Figure 6
gives an overview of all measurements and simulated waveforms of the TDR probe inserted to different depths in water. The simulated waveforms were calculated assuming a two-layer sample. The only parameter that was determined by inverse modeling was the relative length of the probe in the air layer because the Debeye parameters for air and water are known.
View larger version (34K):
[in this window]
[in a new window]
|
Fig. 6. Plots of measured (red circles) and simulated (lines) waveforms of a probe submerged in water to different depths. The top chart is of the 251-point waveform with optimized cable-tester resolution; the bottom chart shows the 1024-point measurements with fixed cable-tester resolution. The waveforms marked with the letter a are difficult to analyze with the algorithm by Heimovaara and Bouten (1990). The first reflection of waveform marked with the letter t is not identified correctly by the tangent analysis method.
|
|
Figure 7
presents the estimated depth in water from the TDR waveforms compared with the depths measured with the high-precision workbench. The insertion depth was measured for each of the three wires of the probe because all three wires had different lengths, with the inner wire being the shortest. The means of the inner and outer wires closely follow the 1:1 line with a positive deviation at small depths and a negative deviation at larger depths. The deviation at the small depths can be explained by the fact that the inner wire is about 7 mm shorter than the two outer wires and that the sensitivity of the TDR probe is largest close to the inner wire. Why the deviation becomes negative as the depth increases is more difficult to explain. It appears that the sensitivity of the outer wires increases because the mean of the three wires crosses the 1:1 at full submergence of the probe (i.e., the estimated depth is exactly the same as the mean length of the wires).
View larger version (9K):
[in this window]
[in a new window]
|
Fig. 7. Plot of the optimized insertion depth in water against the measured mean insertion depth of the mean length of the three wires for Exp. 1.
|
|
Comparison of Results from Inverse Modeling with Direct Travel Time Analysis
In the standard application of TDR for soil water content measurement, the apparent dielectric permittivity of the soil is estimated from the travel time of the TDR pulse in the soil. This travel time is estimated by picking points on the TDR waveform using a tangent method (Heimovaara and Bouten, 1990). Huisman et al. (2002) showed that the analysis of waveforms using the S11 scatter function and inverse modeling gives more accurate results than the direct travel time analysis. We compare the results from the algorithm of Heimovaara and Bouten (1990) with the travel time in each section of our transmission line calculated with:
 | [6] |
at 1 GHz (Lin, 2003b). Figure 8
shows the results of this comparison for the layered water and air system as well as the layered undisturbed soil cores. The travel times cover the complete range from air to water. The points that deviate from the regression line correspond with TDR waveforms in which the time value for the end point as identified with the travel-time algorithm is too large. All four points that fall below the 1:1 line are measurements in a layered system of air and water. The corresponding waveforms are marked in Fig. 6. The errors in the measurements are primarily due to errors in the travel-time analysis with the tangent method. These errors arise because of ambiguous end-point detection for the three waveforms, previously marked with the letter a in Fig. 6. The primary and secondary reflections from the end of the probe are very close to each other, and the algorithm detects a secondary reflection instead of a primary. The other error, marked with the letter t, was caused by a wrong detection of the first reflection. This was caused by noise in the waveform. These errors are quite common for the traditional automated travel-time analysis algorithms.
Statistical evaluation of the parameters of the straight line through these points demonstrates that the slope is not significantly different from 1 on a 95% confidence level, whereas the intercept is significantly different from 0. The proposed inverse modeling method gives higher travel time estimates, thereby resulting in higher values of the apparent permittivities. The intersection point calculated by the algorithm of Heimovaara and Bouten (1990) is primarily based on pragmatism. The calibration for the time domain method apparently causes a slight bias.
The inverse modeling with the SCEM-UA algorithm explicitly accounts for the primary and secondary reflections. As a direct consequence, the accuracy of the SCEM-UA derived results is higher. This is further illustrated in Fig. 9
, which presents a plot of the standard deviation of the optimized water depth, derived from the posterior parameter distribution. The standard deviation is a function of the water depth, and the value ranges from 0.9 µm in air to 29 µm in water. The decrease in accuracy with increasing water depth is mainly caused by dielectric relaxation and conductive losses as the TDR waveform passes through the water sample. The primary and secondary reflections in the waveform become less clear, and as a result, the uncertainty in the inverse modeling approach increases. Another cause for a decrease in accuracy is the fact that the resolution in the measurement of the 251 point waveforms decreases with increasing dielectric permittivity because the software changes the horizontal settings of the cable tester to fit the complete first reflection in 251 points.
View larger version (9K):
[in this window]
[in a new window]
|
Fig. 9. Plot of the standard deviations calculated from the distributions of the optimized insertion depth in water obtained with the SCEM-UA algorithm for Exp. 1.
|
|
Experiment 2: Measurement of the Permittivity Distribution in Packed Layered Samples
Figure 10
shows the first reflections of waveforms measured on the two packed samples. It is clear that the order of packing has a significant effect on the waveforms which may complicate the traditional time domain analysis as it is not immediately clear which reflection in the signal corresponds with the end of the probe. The goal of the inverse modeling of TDR waveforms in this case is to obtain the distribution of the frequency-dependent complex dielectric permittivity along the length of the probe and to use this distribution to assess the water content distribution. We used a sequential approach in which we assumed the sample to consist of two, four, and eight layers. Although we packed the samples in two layers, the high difference in water content between the two layers may cause an immediate redistribution of water in the sample after packing, causing the boundary between the two layers to be diffuse. In addition, the samples were hand-packed and we may expect heterogeneity in density along the probe. Allowing the model to consist of more layers is a way to cope with the diffuse boundary between the layers and the heterogeneities along the sample.
View larger version (13K):
[in this window]
[in a new window]
|
Fig. 10. Plots of the measured waveforms in the two packed samples. The order of packing is marked in the plot and is from probe head to probe end.
|
|
The SCEM-UA derived mean and standard deviations of the parameters for the two-layer optimization on the measurements on both packed layered samples are given in Table 6. Although not further demonstrated here, all marginal posterior parameter distributions were approximately Gaussian. The means of the optimal parameter distributions obtained from the four TDR waveforms cover a wide range: s ranges from 46 to 289, H from 1.8 to 30, frel from 493 Hz (logfrel = 5.693) to 47.5 MHz (logfrel = 7.676), and DC ranges from 0.015 to 0.049 S m–1. The Debeye parameters allow the model to describe dielectric relaxation and as such the parameter values are highly dependent on the frequency bandwidth of the measurement. For our TDR system the bandwidth is estimated to range from 7 MHz to 1.5 GHz [10log(f) from 6.85 to 9.18] (Heimovaara et al., 1996). The very low value for the relaxation frequency for the sand layers is explained by the fact that the sand is a nonrelaxing medium in this frequency range. A very low relaxation frequency causes the dielectric permittivity to be a constant across the frequency band of the measurement. This is a strong indication that no information on relaxation can be found in this specific waveform measurement. The same can be said for very high relaxation frequencies. These do not occur in our results because of our choice of prior ranges. The very high static dielectric permittivity values for the sand layers are also explained by the very low relaxation frequencies. The optimized values are more or less trivial, which is also indicated by the relatively high values of the standard deviation. The standard deviation is a measure of the width of the marginal parameter distribution and is therefore an indication of the information content of the measurement for a certain parameter (Vrugt et al., 2001, 2002; Weerts et al., 2001).
View this table:
[in this window]
[in a new window]
|
Table 6. Mean values for the optimized parameters for the waveform measurements on the packed layered samples from Exp. 2 assuming a two-layer sample and optimizing the length fraction of Transmission-Line Section 1. Standard deviations are given in italics. Layering is indicated in the header.
|
|
Figure 11
shows the comparison between the measured waveforms and the fits after optimization with two, four, and eight layers. Table 7 gives an overview of the sums of squared errors of the different optimizations. It is clear that the quality of the fit improves as more layers are allowed in the model. The waveforms provide enough information to identify the complex dielectric permittivity for models consisting of up to eight layers. Clearly, considering the samples to consist of only two layers is not accurate. The fit could still improve if more layers are included, but the number of iterations required for convergence could become very large.
View this table:
[in this window]
[in a new window]
|
Table 7. Mean sum-of-squared residuals after optimization using the two-, four-, and eight-layer model assumption.
|
|
Figure 12
gives the ranges that the real part of the dielectric permittivity values cover between 100 MHz and 2 GHz as calculated from the optimized parameters for each model layer with Eq. [3]. Thicker lines indicate the occurrence of a relaxation between 100 MHz and 2 GHz. From the figures it is clear the relaxation only occurs in the loam layers and not in the sand layers. This is consistent with the very low relaxation frequency of both sand layers. Wetter, finer-textured soils show more relaxation than drier, coarser soils. Comparison of the top and bottom plots in Fig. 12 shows that similar properties are obtained from the waveform for the sand and loam layers in the different samples. The water content distribution in the sand layer of the loam–sand sample seems to be more heterogeneous that in the sand–loam sample.
View larger version (12K):
[in this window]
[in a new window]
|
Fig. 12. Range of the real part of the dielectric permittivity between 100 MHz and 2 GHz along the probe, calculated from the optimized parameters for the two packed samples. The probe is divided into two sections in the two-layer model. For the four- and eight-layer models the number of sections was doubled by halving the lengths of the two- and four-layer models.
|
|
The frequency-dependent complex dielectric permittivity is described with Eq. [3] and consists of four parameters, including DC. Dielectric relaxation and other losses along the sample reduce the magnitude of the higher frequencies in the voltage wave passing through the sample. As a result the effective frequency bandwidth changes as the voltage wave moves through the sample. Traditional TDR measurements of soil water content are based on Ka, which is calculated from the travel time in the probe. As we do not exactly know which frequency should be used to assess the apparent dielectric permittivity, we chose the real part of the complex dielectric permittivity at 1 GHz as an estimate of Ka (Lin, 2003b). In doing this we may underestimate the permittivity for situations where the frequency bandwidth does not contain these high frequencies.
Figure 13
gives an overview of the 90% confidence interval of the estimated water content using Eq. [5] for the three models. The 90% confidence interval is estimated by sampling the optimized parameter distribution. The results indicate that the assumption of a smoothing of the water content gradient in the sample is reasonable for the loam–sand sample. The gradient in the water content is still rather steep because the two- and four-layer models are very similar.
View larger version (11K):
[in this window]
[in a new window]
|
Fig. 13. Ninety-percent confidence interval for the water content distribution in the packed samples for the two-, four-, and eight-layer models. Probability was calculated by sampling the optimized parameter distributions.
|
|
Figures 12 and 13 provide us with insight in the quality of this inverse modeling procedure. It is important to realize that in this approach we assume that the models we use are accurate and that any model errors can be partly compensated for in the optimized model parameters. Given the quality of the simulated wave forms in water and air after calibration, we believe that model errors in the description of our measurement setup (cable tester, cables, and probe) are minor. This was also proven by the results from the experiment in water and air. In addition to the uncertainty in the description of the measurement setup, the uncertainty in water content is also related to the models used to describe the dielectric permittivity in the sample and the translation from the dielectric permittivity to the water content. We assume that the complex dielectric permittivity of our sample can be described with a simple Debeye relaxation given by Eq. [3]. This model allows for a single relaxation frequency in the frequency bandwidth of measurement. This model can be too simple for soils as it cannot account for multiple relaxation processes within the bandwidth such as the Maxwell–Wagner relaxation due to double-layer or bound-water relaxation (Hilhorst, 1998). Topp's equation, Eq. [5], is used to calculate the water content. To apply this equation we estimate the apparent permittivity with the real part of the complex permittivity at 1 GHz. Figure 12 shows that for certain layers the permittivity seen in the frequency bandwidth can show a significant variation. Using a wrong frequency in Eq. [6] can result in significant error. In addition to these uncertainties due to model errors, we also have the parameter uncertainty in the optimized parameters. This uncertainty found in the parameter distribution is rather small and is illustrated by the width of the lines in Fig. 13. The uncertainty due to the unknown frequency, and an error in the model for the dielectric properties is not accounted for in this bandwidth. As a result the true uncertainty is larger.
Table 8 gives an overview of the water content values calculated from the gravimetric water content of the loam and sand and thickness of the packed layers and the water content ranges obtained from the inverse modeling. The estimated values are close to the values measured gravimetrically, especially considering that no soil specific calibration has been used and that water redistribution can have taken place in the samples, causing the water content of the sand layer to increase.
View this table:
[in this window]
[in a new window]
|
Table 8. Water content and dry bulk density values of the layers in the packed samples calculated from layer thickness (dL) and the gravimetric water content of the loam and sand. The gravimetric values were determined by oven drying after the experiments were conducted.
|
|
|
SUMMARY AND CONCLUSIONS
|
|---|
We presented a combined MSSF model and stochastic global optimization method (SCEM-UA) for interpreting measurements from an arbitrary (field-type) TDR setup. The method enables the identification of water saturation profiles along a TDR probe. The approach can be applied to standard TDR probes, and no special input waveforms need to be measured. As such, the waveform interpretation is completely automated, which increases the reproducibility and accuracy of the TDR method because all information in the TDR waveform from primary and secondary reflections in the waveform is used in this approach.
An accurate calibration of the TDR measurement system is essential because a misfit in the calibration can result in erroneous parameters obtained by inverse modeling. We found that a three-step calibration provided good results. A waveform of a 50- load was used to identify the characteristics of the input signal. Measurements with open-end cables were used to obtain an initial estimation of the parameters used to describe the internal circuitry of the cable tester, the connectors, and cables. Finally, using the characteristics of the input signal and the initial estimation of parameters, the measurement setup was calibrated on measurements in air and water with a distinct electrical conductivity. Knowledge of the electrical conductivity of the water is not required.
The accuracy of the method is illustrated with the results of a lab experiment in a layered system of water and air. The standard deviation of the distribution of the optimized values of measured water depths ranged from 0.9 to 29 µm. The value of the standard deviation increases with the water depth and can be explained by the fact the reflections in the TDR signal become less clear due to dielectric and conductive losses. This result is found for all measurements performed in our study. A consequence of this result is that the accuracy of the TDR measurement is limited by probe length, a result already well known (Robinson et al., 2003). Theoretically the method is applicable to probes with arbitrary lengths, but experimental work is required to test the maximal feasible length.
Measured waveforms on layered soil samples were analyzed with the SCEM-UA inverse modeling scheme using models for the TDR probe that consisted of one to eight layers. The sum-of-squared residuals decreased with increasing number of layers, but the number of iterations required for convergence increased as well. The approach with eight layers clearly was able to accurately describe measured waveforms. The water content distribution obtained from the results showed a large variation in water content within the different layers of the sample. The smallest layer for which we demonstrated the possibility to determine the dielectric properties and thus the associated water content was approximately 0.012 m.
We showed that we can accurately describe the TDR measurement setup. We assumed a rather parsimonious single-relaxation model. Our results indicate that using the single Debeye relaxation given in Eq. [3] is probably too simple for describing the dielectric properties of (wet) soils. Recently, however, other mixing models have become available, which are based on the dielectric properties of the soil constituents and the physical structure of the porous medium (i.e., Lin, 2003b; Robinson and Friedman, 2003). In addition, a large body of theoretical work has accumulated on effects that air gaps have on the TDR measurement. These models and approaches can easily be built into the methodology we presented. If an accurate model is found, TDR might also prove to be useful for determining other properties of the soil such as the dry bulk density, which is an important property in volumetric mixing models.
The availability of the proposed inverse modeling framework increases the applicability of the TDR method. Because of its high accuracy and potential for fully automated data analysis, the method is suited to obtain detailed (real-time) data on heterogeneous infiltration.
|
APPENDIX
|
|---|
The Multi-Section Scatter Function Model
Feng et al. (1999) presented a multisection S11 (reflection) scatter function to model the waveform measured by a TDR cable tester. In this paper we adopt the same numbering convention as proposed by Feng et al. (1999), as also shown in Fig. 2. The corresponding multisection S11 scatter function is given by Feng et al. (1999) as
 | [A1] |
with
 | [A2] |
being the reflection coefficient between the different sections k of the transmission line indicated in Fig. 2, Zk(f) the impedance and 
k(f) the propagation coefficient of section k and Ll,k the length of each section k. Equation [A1] is calculated by iterating through all sections, starting with Section 1. To calculate the scatter function of a multisection transmission line we must make an assumption concerning the final reflection (i.e., the reflection at Transmission-Line Section 0). For an open-ended transmission line, S110 = 1; for a shorted transmission line, S110 = –1; and for a matched transmission line S110 = 0. The characteristic impedance of the last section in our transmission line (i.e., k + 1 = 7) is assumed to be 50 .
The parameters Zk(f) and k(f) can be calculated from models derived from transmission-line theory. Models are available for nonideal coaxial transmission lines in which we can account for the distributed resistance Rs along the conductors, the inductance L, the shunt conductance G, and the capacitance C per unit length of the transmission line. It is possible to calculate Z(f) and (f) from the geometry of each coaxial transmission-line section using the following equations (Ram
TOP
ABSTRACT
INTRODUCTION
