Published online 25 February 2008
Published in Vadose Zone J 7:238-248 (2008)
DOI: 10.2136/vzj2007.0087
© 2008 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SPECIAL SECTION: GROUND PENETRATING RADAR IN HYDROGEOPHYSICS
Identifying Unsaturated Hydraulic Parameters Using an Integrated Data Fusion Approach on Cross-Borehole Geophysical Data
Majken C. Loomsa,*,
Andrew Binleyb,
Karsten H. Jensena,
Lars Nielsena and
Thomas M. Hansenc
a Univ. of Copenhagen, Dep. of Geography and Geology, Øster Voldgade 10, DK-1350 Copenhagen K, Denmark
b Lancaster Univ., Dep. of Environmental Science, Lancaster, LA1 4YQ, UK
c Univ. of Copenhagen, Niels Bohr Institute, Juliane Maries Vej 28, DK-2100 Copenhagen Ø, Denmark
* Corresponding author (mcl{at}geol.ku.dk).
All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.
Received 11 May 2007.
 |
ABSTRACT
|
|---|
Cross-borehole geophysical data can provide valuable information concerning hydrologic properties of the unsaturated zone. Such data are most often used sequentially, where images of soil physical properties are obtained through numerical inversion and then converted to hydrologic state properties using petrophysical relationships. If not accounted for, inversion artifacts are transferred to the resulting hydrologic images. We propose a framework in which multiple geophysical data sets can be incorporated using an integrated data fusion approach. The geophysical data collected are integrated in a forward modeling approach to evaluate a series of plausible hydrologic models. The approach permits an evaluation of the sensitivity of geophysical data for constraining hydrologic model parameters. We illustrate the approach using geophysical data collected during a dual water and solute tracer experiment in the unsaturated zone. Cross-borehole ground penetrating radar and electrical resistivity tomography, measuring electromagnetic travel time and electrical transfer resistances, respectively, were collected during a 20-d period. As a first approximation, one-dimensional flow was considered and three models (one, two, and five layers) of the subsurface were evaluated. The five-layered model was found to be the only model capable of mimicking the infiltration pattern satisfactorily. The results showed that only the hydraulic conductivity and one of the parameters (empirical parameter n) describing the soil moisture release curve for three of the five layers could be constrained by the data, illustrating the nonuniqueness of the problem. The data fusion approach proved, however, that application of multiple geophysical methods may reduce hydraulic parameter uncertainty.
Abbreviations: ERT, electrical resistivity tomography • GLUE, Generalized Likelihood Uncertainty Estimation • GPR, ground penetrating radar • ZOP, zero-offset profile
|
INTRODUCTION
|
|---|
Cross-borehole geophysical methods are increasingly used to estimate hydrogeologic properties of the unsaturated zone. These methods provide high-resolution in situ information for large support volumes of the subsurface with minimal intrusion and thereby have obvious advantages over other more traditional methods such as neutron probe and time domain reflectrometry.
In particular, cross-borehole electrical resistance tomography (ERT) and ground penetrating radar (GPR) have become very popular in the past decades. Applications of these methods range from characterization of permeable pathways in fractured basalt (Hubbard et al., 1997) and monitoring preferential flow during early stages of snowmelt (French et al., 2002) to following the effect of different remediation processes (Daily et al., 1995; LaBrecque et al., 1996). The moisture content estimates derived from GPR have been shown to have a high reproducibility and accuracy (Parkin et al., 2000; Ferré et al., 2003). As a result, water front movement caused by natural infiltration can be monitored successfully despite just slight variations in moisture content (Binley et al., 2002b).
Cross-borehole GPR and ERT are especially suitable for monitoring relative changes in soil moisture and electrical conductivity, since local variations in the petrophysical relationships are, in this way, reduced or completely avoided (Binley et al., 2002a). Therefore many studies have been performed using the two methods to investigate the flow dynamics of the unsaturated zone resulting from forced loading to the system (e.g., Daily et al., 1992; Binley et al., 2001; Alumbaugh et al., 2002). Water flow patterns have been visualized and attempts have also been made to quantify important flow and transport parameters, especially the effective hydraulic conductivity (Ks) (Binley et al., 2002a; Cassiani et al., 2004: Kowalsky et al., 2004, 2005; Rucker and Ferré, 2004; Cassiani and Binley, 2005; Cassiani et al., 2006). Estimation of subsurface hydraulic parameters can be performed by combining moisture content data from cross-borehole GPR with collected bulk electrical resistivity data obtained from cross-borehole ERT data (Winship et al., 2006). Looms et al. (2008) applied the two methods simultaneously for determination of fluid conductivity and thereby tracer concentration, which was subsequently used for moment analysis of the evolving plume.
The tomographic images resulting from cross-borehole geophysical surveys are, however, subject to several uncertainties and artifacts caused by the inversion techniques used to invert the collected data. The inverse problem often consists of a mixed-determined problem, and to solve this, regularization (e.g., Constable et al., 1987) of the problem is normally performed. As a result, the images are typically smooth, minimum-variance estimates of the subsurface parameters and cannot capture the fluctuations of the true model parameters (Hansen et al., 2008). Varying resolution of the model parameters can also affect the images and the results (Singha and Gorelick, 2005), and model resolution is an important factor to consider if the inversion images are used for further analysis. In addition, combining geophysical images from multiple methods in any analysis is nontrivial because of the method-specific resolution and sensitivity, as illustrated by Day-Lewis et al. (2005).
Various techniques have been proposed to account for the varying resolution of the tomographic images, including random field averaging (Day-Lewis et al., 2005), apparent petrophysical relations (Singha and Gorelick, 2006; Singha and Moysey, 2006), and full-inverse statistical calibration (Moysey et al., 2006). Also, joint inversion of different types of geophysical data with complementary resolution properties can help produce more reliable geophysical models (Kowalsky et al., 2005; Linde et al., 2006).
Using an "integrated data fusion" approach, the geophysical data are not inverted to create images of the subsurface but are solely used to evaluate the likelihood of a series of plausible hydrologic models with varying unsaturated hydraulic parameters. This approach was used with one-dimensional GPR data collected during natural infiltration in previous work by Binley and Beven (2003). The Generalized Likelihood Uncertainty Estimation (GLUE) method (Beven and Binley, 1992) was used as a basis for this approach. They found that it was difficult to determine the hydraulic parameters based on the data collected, as the natural changes in moisture content during the data collection period were very small. Even multiple profiles with time appeared to offer very little additional information.
In this study, we applied the integrated data fusion approach presented previously to a new data set describing large dynamic changes in unsaturated sand. The data were collected during a forced infiltration experiment using both solute and water as tracers. To resolve the concentration changes and the moisture content variations, two different data types (cross-borehole ERT and GPR) were included in the analysis. The sensitivity of the two data types to constrain the unsaturated hydraulic parameters of the subsurface and the limitations connected to the data fusion approach can thereby be evaluated.
|
Field Site and Field Experiment
|
|---|
A field site was established at Arrenæs, Denmark, on a 20- to 30-m unsaturated deposit of alluvial sediments, consisting mainly of sands. In June 2004, the field site was instrumented with eight 12-m-deep boreholes with the purpose of studying the flow and transport processes through the unsaturated zone. The eight boreholes formed a cross consisting of two lines (see Fig. 1
). Along each line, the two outer boreholes (7 m apart) were equipped with ERT-instrumented polyvinyl chloride (PVC) tubes (electrodes every 50 cm), while the two inner boreholes (5 m apart) had PVC access tubes for GPR antennae.
View larger version (29K):
[in this window]
[in a new window]
|
FIG. 1. Schematic drawing of the field site setup at Arrenæs, Denmark. The light gray area indicates the infiltration area.
|
|
Initial GPR and ERT measurements showed a very low seasonal variation in moisture content at the field site. A forced infiltration experiment was therefore conducted in the autumn of 2005 to create dynamic flow conditions in the subsurface. During a period of 20 d, 
95,000 L of clean water was irrigated through 484 drippers over a 7.33- by 7.33-m area. After 4 d of infiltration, 890 L of saline tracer was added through the drippers. The saline tracer had an electrical conductivity of 10.54 S m–1 and was produced by adding 75 kg NaCl to 890 L clean water (76 mS m–1).
During the field experiment, 14 measurement surveys were performed, including a background (pretracer and preinfiltration) survey. The background moisture content profiles estimated by GPR and ERT, shown in Fig. 2a
, were found to be consistent, exhibiting similar trends throughout the measurement depth. Furthermore, the preinfiltration moisture content profiles were found to be in good agreement with grain size analysis results from samples extracted from a nearby borehole, Fig. 2b, since high values of moisture content were observed at depths containing fine sediments and vice versa.
View larger version (22K):
[in this window]
[in a new window]
|
FIG. 2. (a) Background moisture content profiles estimated using cross-borehole ground penetrating radar (GPR) and electrical resistivity tomography (ERT); (b) sediment samples from a nearby well (d10, d50 and d90 are the 10th, 50th and 90th percentiles of the grain size distribution; data provided by Copenhagen Energy), and (c) the five-layered model used in the hydrologic forward simulation. The boundaries of the different materials were deduced from (a) and (b).
|
|
Moisture content and tracer mass profiles estimated from inverted cross-borehole GPR and ERT data collected during the 20-d experiment showed that large quantities of water were diverted out of the area. During the first 8 d, a 50% water loss was observed. The consistency of this observed trend has led us to believe that the observed water loss was caused by lateral flow occurring just below the ground surface, and the collected data therefore describes flow patterns resulting from half the applied irrigation. For more details of the infiltration experiment and results, see Looms et al. (2008).
|
Methodology
|
|---|
Richards' equation is often used to describe variably saturated flow in porous media. One-dimensional vertical flow can be expressed as
 | [1] |
where z is elevation, t is time, 
(h) is the volumetric moisture content, and K(h) is the hydraulic conductivity. We assume here that the retention and hydraulic conductivity functions can be represented by the parametric models of van Genuchten (1980):
 | [2] |
 | [3] |
 | [4] |
 | [5] |
where r is the residual moisture content, s is the saturated moisture content (or porosity), 
and n are empirical parameters, and Ks is the saturated hydraulic conductivity.
A total of five model parameters (r, s, , n, and Ks) therefore describe the hydraulic properties for each soil layer considered. We used the methodology depicted in Fig. 3
to examine the performance of different sets of model parameters by comparing forward calculated responses of geophysical measurements, corresponding to the specific flow conditions determined by the model parameters, with the geophysical measurements collected in the field. Below, each step of the process is described individually.
Solute transport was described by convective–dispersive transport (
im
nek et al., 1998). We assumed fully conservative transport with a longitudinal dispersivity of 25 cm, representing approximately 1/10 the length scale of the tracer migration during the first 4 d. An average dispersivity value, estimated using moment analysis on the inverted geophysical data, supports this value (Looms et al., 2008).
Sampling Parameter Values
The soil profile was discretized into a number of horizontal layers (the details of which are described below). For each layer, the hydrological parameters r, s, , n, and Ks need to be defined to simulate water migration through the profile. Table 1
lists parameter ranges determined from retention experiments on core samples taken from the field site, as well as published ranges for sand from the ROSETTA database (Schaap et al., 2001). Table 1 also lists the parameter ranges (a priori ranges) selected for the simulations. Note that the parameter ranges of , n, and Ks have been selected as widely as possible in order to not exclude plausible parameter values, and the minimum value of hydraulic conductivity has been lowered several orders of magnitude to better describe potential fine-grained materials near the surface. Furthermore, ERT and GPR moisture content profiles indicated background moisture contents down to 0.07 (see Fig. 2a). As a result, the range of the residual moisture content, r, was limited to <0.06 to avoid negative effective saturations. Finally, as porosity exhibited a low degree of variation, this parameter was fixed at 0.41 to limit the number of varying parameters. The selected value represents the average measured porosity of sand at 1- to 2-m depth.
View this table:
[in this window]
[in a new window]
|
TABLE 1. Range of van Genuchten parameters estimated from core samples collected at the field site; published for sand found in the ROSETTA database, computed from the average and one standard deviation of 308 sands; and selected for the HYDRUS-1D simulations.
|
|
In this study, the hydraulic parameters were assumed to be uncorrelated. In future work, it would be of benefit to include correlation between parameters to constrain the parameter ranges, thus allowing more efficient sampling of the parameter space. Here, the parameter ranges were sampled using a Latin hypercube method. This method works like a constrained Monte Carlo sampling approach (McKay et al., 1979). Each parameter range is subdivided into M equally probable intervals. For each interval, one value is selected at random and subsequently paired at random with similarly determined values of the remaining parameters, resulting in M independent realizations. All parameters were assumed to be uniformly distributed within the specified intervals. The hydraulic conductivity, Ks, was sampled in the log-transformed space. A public domain software Latin hypercube sampling tool was used (www.mathepi.com/epitools; verified 19 Dec. 2007).
Hydrologic Model
In the next step, the selected hydrologic parameters were used as input to the one-dimensional unsaturated flow and transport code HYDRUS-1D (imnek et al., 1998). Previous analysis of the geophysical data (Looms et al., 2008) indicated that a one-dimensional assumption is probably not true in the shallow horizons, as significant lateral flow is likely to have occurred as a result of an impeding layer near the soil surface. We, therefore, ignored this lateral flow and adopted a reduced "effective" infiltration rate to maintain the total tracer mass observed in the deeper soils by the geophysical data.
Flow and transport simulations were run for each parameter set (realization) for the same initial and boundary conditions. The upper boundary condition was a specified flux of 0.1842 cm h–1, corresponding to half of the applied irrigation. After 4 d of infiltration, the saline tracer was applied over 9 h. The lower boundary condition of the model was in the form of atmospheric pressure at 30-m depth corresponding to the known position of the water table. The vertical discretization of the model was 1.25 cm throughout the profile, resulting in a total of 2401 node points. The average background moisture content profile collected using the cross-borehole GPR and ERT methods (Fig. 2a) was used as the initial condition for the model. Evapotranspiration was not considered since calculations indicated that the evapotranspiration would constitute only an insignificant fraction of the applied irrigation (<1%).
Moisture content and tracer concentration profiles were computed for eight selected times (Days 0, 1, 2, 3, 4, 5, 6, and 8) matching the first eight ERT and GPR measurement times.
Geophysical Forward Analysis
The eight moisture content and concentration profiles calculated by HYDRUS-1D were subsequently converted to electromagnetic velocity (v) and bulk electrical resistivity (
b) profiles using well-established petrophysical relationships. Site-specific parameterization of these relationships is reported in Looms et al. (2008). The resulting moisture content and concentration profiles were then used in forward numerical solvers to calculate electromagnetic travel times and electrical transfer resistances corresponding to the same measurement configurations as applied in the field.
Constant petrophysical parameters were assumed throughout the investigated soil column, since the sediments in the soil volume consist mainly of sands of varying coarseness, with clay contents <1% and negligible organic contents. The topsoil (approximately 1.0–1.5 m thick) does, however, contain high clay contents (14%) and a bias may therefore be present in the ERT-estimated resistivities representing this layer. The GPR measurements collected within the top 1.5 m were not used in the analysis, as will be explained below.
Ground Penetrating Radar
In low-loss material and at high frequency, the electromagnetic velocity can be estimated from (e.g., Reynolds, 1997)
 | [6] |
where 
r is the relative permittivity of the porous medium and c is the radar wave velocity in air (0.3 m ns–1).
The relative permittivity is related to moisture content and we used the empirical formula of Topp et al. (1980):
 | [7] |
The velocity distribution was subsequently used to compute GPR travel times for a given set of measurements. For this purpose, the nonlinear partial differential eikonal equation was solved using the FAST code (Zelt and Barton, 1998). A bending-ray algorithm was chosen instead of a straight-ray algorithm due to the potentially large velocity contrasts arising at the wetting front.
Electrical Resistivity Tomography
The bulk electrical resistivity was calculated using Archie's law (Archie, 1942):
 | [8] |
where w is the pore water electrical resistivity, 
is the porosity, and m and n are empirical constants typically set to 1.3 and 2, respectively, but dependent on the soil in question.
The three-dimensional finite-element-based electrical resistivity code R3 (Binley, 2007) was used to calculate electrical transfer resistances corresponding to the actual measurements collected. The region was discretized into 59 by 52 by 52 cells having a vertical, horizontal, and transverse discretization of 0.25 m in the infiltration area. This setup resulted in numerical errors that were typically around 1%.
Geophysical Field Data
Each field measurement survey consisted of 1767 ERT and 43 GPR measurements. The ERT data were collected using 96 electrodes: 84 borehole electrodes installed in four ERT boreholes and 12 surface electrodes. Only horizontal borehole dipoles were used, resulting in higher values of measured transfer resistances (Winship et al., 2006).
The GPR data were collected using the zero-offset profile (ZOP) mode. Two antennae were lowered simultaneously into a set of boreholes, stopping every 0.25 m to obtain a measurement. The data collected from 0.00 to 1.50 m were later disregarded. In this region, the first arriving electromagnetic signal was a refracted wave travelling at the air–soil interface and not the direct wave through the ground. A ZOP therefore consisted of 43 measurements.
For a more extensive description of the measurements conducted and the resulting inversion results using traditional algorithms, see Looms et al. (2008).
Measurement of Data Misfit
To evaluate the extent to which the individual hydrologic models represent field conditions, the measured and simulated data were compared using a misfit value:
 | [9] |
and
 | [10] |
where r and tt are the average relative data error of the transfer resistances (=0.08) and the absolute data error of the travel times (=0.8 ns), respectively. The data errors represent various error sources, e.g., measurement errors, modeling errors, and errors associated with picking the first arrival of the direct electromagnetic wave. The values used in this analysis are consistent with data error estimates used for traditional inversion techniques (see Looms et al., 2008).
As described above, the HYDRUS-1D model was initiated using the average ERT and GPR background moisture content profile. The electrical resistance and travel time measurements collected at Day 0 were therefore not identical to the simulation values of Day 0. As a result, a misfit value existed for Day 0 (i.e., MISFITERT = 1.75 and MISFITGPR = 2.07). These misfit values are the minimum misfits expected between the simulated and measured data. The measurement types (ERT and GPR) were therefore normalized according to these background misfit values in the combined analysis. When data from multiple days were used, a MISFIT value was calculated for each day (according to Eq. [9] and [10]), after which each day was equally weighted. It should be noted that these were all subjective choices that may have influenced the results of the performed analyses. Additional misfit values could be estimated using different weights and misfit functions; however, an extensive analysis into this topic is beyond the scope of this work.
|
Results
|
|---|
A homogeneous profile was considered first, assuming uniform hydraulic parameters throughout the investigated domain, i.e., a one-layered model; however, close inspection of grain size analyses from core sampling at a nearby borehole, as well as one-dimensional profiles of background soil moisture contents (Fig. 2), indicated that the subsurface may be subdivided into two or even five soil layers with different hydraulic characteristics. Two additional systems were therefore also tested: a two-layered system, where the profile is divided into an upper 1.5-m horizon and an underlying uniform horizon, and the five-layered setup illustrated in Fig. 2c.
Initially, the proposed methodology was used to assess the different hydrologic realizations using only GPR data because of the significantly lower computational cost of the GPR forward modeling compared with that of the ERT problem. Furthermore, all eight measurement days were used in this initial analysis. To reduce the computational costs of the problem, however, we only included data from Day 8 in the combined ERT and GPR data analysis.
To evaluate whether the unsaturated hydraulic parameters are constrained by the data, we plotted the misfit value (calculated from Eq. [9] and [10]) as a function of parameter values (r, , n, and Ks) for the realizations tested (Fig. 4
, 6, 9, 10, and 11). For illustration purposes, only the 500 realizations with the lowest misfit are included in the figures. Furthermore, we highlighted the range of the five best misfit results with a red line. By comparing the extent of this range with the a priori range of the selected parameter and by evaluating whether the dots with low misfit values are clustered at certain parameter ranges, it becomes clear whether the data can reduce the uncertainty of the hydraulic parameters investigated.
View larger version (36K):
[in this window]
[in a new window]
|
FIG. 9. A synthetic evaluation of the five-layered model analysis using the two data types combined. For illustration purposes, only the 500 realizations with the lowest misfit values are included in the figure. The parameter ranges of the five realizations with the lowest misfit values are indicated with red and the "true" parameter value is shown with a red diamond.
|
|
View larger version (31K):
[in this window]
[in a new window]
|
FIG. 11. The misfit values of the five-layered model using only electrical resistivity tomography (ERT) data, using only ground penetrating radar (GPR) data, and combining ERT and GPR data in the analysis. For illustration purposes, only the 500 realizations with the lowest misfit are included in the figure. Only the misfit plot of the top three layers of the saturated hydraulic conductivity, Ks, and the empirical parameter n are shown.
|
|
One-Layered Model
Figure 4 depicts the computed misfit between the modeled and the measured data (Eq. [10]) according to the parameter values used in 1000 individual one-layered model realizations. The parameter ranges of the five best results, i.e., lowest misfit values, are tabulated in Table 2
. It is evident, from Fig. 4 and Table 2, that the range of the parameter n is not constrained by the data, since the range of the five lowest misfit realizations is almost identical to the a priori range and the dots do not exhibit any structure; however, the parameter ranges of r, , and Ks are somewhat constrained compared with the a priori range. Values of Ks, centered around 0.40 cm h–1, appear to result in slightly lower misfits between the simulated and measured travel times, and this trend is supported by the misfit plot that exhibits a clear clustering at this intermediate Ks value.
To illustrate the variance of the 1000 realizations and how the five realizations with the lowest misfit values compare, the simulated travel time data for all the realizations are plotted in Fig. 5
for three selected days. The ZOP travel times give a one-dimensional representation of the water movement with time, since travel time recordings are directly proportional to the moisture content (Eq. [6–7]). The five best misfit results have been highlighted with red lines and the measured travel time profile is plotted in green. Figure 5a shows the travel time profiles of Day 0, representing the initial condition in the hydrologic model before water was added to the system. The simulated moisture content profile is therefore the same for all realizations. After just 1 d of infiltration, the simulations show considerable diversity (Fig. 5b). The five best simulations are, however, almost identical and resemble the measurements very well. Figure 5c shows the last measurement time used in this analysis, Day 8. In this plot, the five best simulations are no longer consistent with the measured travel times. These simulations were based on low saturated hydraulic conductivities and as a result the water front was migrating slowly through the subsurface, increasing the travel times down to only the 2.50-m depth. The increased moisture content, from 3.00 to 8.00 m measured by the GPR method at Day 8, is thereby not captured by these realizations.
View larger version (32K):
[in this window]
[in a new window]
|
FIG. 5. The simulated zero-offset profile (ZOP) travel times resulting from the one-layered model for three selected days: (a) Day 0, (b) Day 1, and (c) Day 8. The five realizations with the lowest misfit values are highlighted with red. The measured travel time profile is highlighted with green.
|
|
A similar analysis can be made using only the travel time data from Day 8 to calculate the misfit values. The results of this analysis are presented in Fig. 6
and 7
. The low misfit values plotted in Fig. 6 appear to have slightly more structure than observed in the similar Fig. 4, suggesting that the misfit values are more parameter dependent. Low misfit values are observed at low values of r (<0.02) and high values of n (>4.15). Furthermore, low values of Ks result in low misfit values. Apart from one of the five best realizations where the Ks value = 7.17 cm h–1 (indicated with an arrow), the remaining four realizations have Ks values <0.16 cm h–1; however, only the realization with Ks = 7.17 cm h–1 (again indicated with an arrow in Fig. 7b) is seen to result in the water front reaching 8-m depth at Day 8, as was observed in the measured travel time profile. The travel times of this realization do not, however, fit the initial moisture contents at Day 1 very well (Fig. 7a).
View larger version (51K):
[in this window]
[in a new window]
|
FIG. 7. The simulated zero-offset profile (ZOP) travel times resulting from the one-layered model for two selected days: (a) Day 1 and (b) Day 8. Only data collected from Day 8 were used to calculate the misfit values. The five realizations with the lowest misfit values are highlighted with red. The measured travel time profile is highlighted with green. An arrow points out the realization with a higher saturated hydraulic conductivity value.
|
|
The fact that the realizations that fit the travel times at Day 8 do not fit the early measurement times and vice versa leads us to conclude that a simple one-layered model is not appropriate to model the observed water infiltration.
Two-Layered Model
To evaluate whether the low grain size material in the top 1.50 m caused the discrepancies observed above, a two-layered model was tested. Again 1000 realizations with different hydrologic parameter sets were run in HYDRUS-1D. For both layers, the four hydraulic parameters (r, , n, and Ks) were varied independently within the tabulated ranges in Table 1. The results of the analyses (not shown here) did not appear to improve the fit to the measured data significantly compared with the one-layered model.
Five-Layered Model
To capture the dynamics of the infiltration pattern, it was necessary to elaborate the conceptual geologic model further to a five-layered model. Again 1000 realizations were generated, each layer varying the four hydraulic parameters (20 parameters in all) within the ranges shown in Table 1. The results of the ZOP travel time profiles, with the five lowest misfits highlighted in red, are plotted in Fig. 8
. Two of the five most optimal realizations now exhibit much better visual agreement with the observed flow patterns, as well as a lower misfit value (i.e., 2.14).
View larger version (54K):
[in this window]
[in a new window]
|
FIG. 8. Travel time profiles from the five-layered model for two selected days: (a) Day 1 and (b) Day 8. The five realizations with the lowest misfit values are highlighted with red. The measured travel time profile is highlighted with green.
|
|
Ground Penetrating Radar and Electrical Resistance Tomography Data Combined
Due to the improved comparison between simulated and measured GPR travel times, the more complex five-layered model was selected for further analysis using the cross-borehole ERT and GPR data combined. A total of 4000 realizations of the five-layered model were run.
Initially, we used one of the realizations as assumed measured ERT and GPR data values to test the limitations of the proposed methodology. In this way, we were able to evaluate whether the misfit analysis was able to identify the "true" parameters unequivocally. To mimic the problem at hand as closely as possible, we chose the realization with the lowest misfit in the combined analysis, listed in Table 3
, as the synthetic model. The results are presented in Fig. 9
, where the "true" parameter values are indicated with a red diamond.
The dots in Fig. 9 and the range of the five lowest misfit values indicate that the proposed methodology can constrain the hydraulic conductivity of Layers 1, 2, and 3 and give some constraint on the moisture release curve parameter n of Layers 1, 2, and 3. There is an especially high clustering of the dots corresponding to the Ks values of Layer 2. The data do not provide enough information to constrain the remaining parameters. Here all the dots appear randomly distributed across the a priori parameter range and the correct parameter value (red diamond) is not included within the range of the five lowest misfit values for three of the varied parameters. An attempt to determine these parameters through conventional optimum searching approaches would therefore lead to nonunique results. Since the water front of the chosen realization did not penetrate into Layers 4 and 5 (8 m), the data cannot constrain the hydraulic parameters of these layers.
To evaluate the effect of using data sets from multiple days in the analysis, we present misfit plots using only GPR for Day 8 and all 8 d in Fig. 10a and 10b
, respectively. The misfit plots show that the additional data from Days 1 to 6 only result in a narrowing of the estimated Ks ranges of Layers 1 and 2, while the parameter range of Layer 3 is increased. Again, the geophysical data were observed to not contain information regarding the bottom two layers, i.e., Layers 4 and 5.
Since the synthetic modeling exercise highlights that we are only likely to be able to constrain parameter Ks and n and, furthermore, this is limited to Layers 1, 2, and perhaps 3, we restricted the presentation of findings from our field data to these parameters of the model (Fig. 11
). In Fig. 11, we further investigate the two data types (ERT and GPR) separately, as well as in a combined analysis.
The ERT data collected in the field provide a means of somewhat constraining the saturated hydraulic conductivity values of Layers 1, 2, and 3. Clear structures in the computed misfit scatter plots of Layers 1 and 2 indicate that values of Ks < 0.1 cm h–1 are not very likely to represent the hydraulic properties of the top 4 m. The ranges of the parameter n are, however, poorly constrained by the ERT data.
The GPR data also provide information to constrain the Ks values of Layers 1, 2, and 3 (Fig. 11). Low Ks values are not very likely, corroborating the finding using the ERT data. Specifically, the Ks values of Layers 2 and 3 are constrained above 1 cm h–1, whereas the hydraulic conductivity value of Layer 1 appears to be more moderate, around 1 cm h–1. Even though GPR measurements were not performed in Layer 1, the Ks value of this layer is somewhat constrained, although to a lesser degree than in Layer 2. In Layer 2, the misfit plot of the van Genuchten parameter n, furthermore, exhibits some structure, suggesting that values of n < 2 are more improbable than higher values (suggesting a well-sorted sand layer).
The misfit scatter plots were not improved dramatically when data from the two geophysical methods were combined (see Fig. 11). The constraints of the parameters, in particular the Ks values of Layers 1 and 2, however, are improved slightly, suggesting that the hydraulic conductivity values of Layers 2 and 3 are higher than that of Layer 1. This was expected, since the estimated grain size distributions for different depths (Fig. 2b) indicated coarse sand below 1.5 m of fine topsoil.
The travel time profiles of the five best realizations of the combined analysis are shown in Fig. 12
. Even though the profiles mimic the overall behavior of the measured travel time profiles satisfactorily, a closer inspection of the simulation profiles reveals secondary water fronts moving through the subsurface during the initial measurement days (indicated with arrows in Fig. 12). The discrepancies become especially evident when the five moisture content profiles are presented (see Fig. 13
). The local spikes in moisture content observed in Fig. 13 were caused by water draining from the initial moisture content profile, shown in black, and therefore represent the inadequacy of the relatively simple five-layered model to conceptualize the heterogeneous subsurface. The small-scale moisture content changes observed in Fig. 13 are, in the velocity profiles (Fig. 12), smoothed out by limited vertical resolution or bending rays travelling preferentially in high-velocity areas.
View larger version (26K):
[in this window]
[in a new window]
|
FIG. 12. Travel time profiles as they developed with time. The five best realizations of the combined analysis (red) are shown along with the measured profiles (green).
|
|
View larger version (38K):
[in this window]
[in a new window]
|
FIG. 13. Moisture content profiles of the five best realizations of the combined analysis as they developed with time. The background moisture content profile (Day 0) is illustrated with black in the subsequent seven measurement times, while the range of all the realizations of moisture content curves is indicated with gray.
|
|
Finally, the five best relative concentration profiles are presented (see Fig. 14
). Four of the profiles are nearly identical, and only one realization is slightly shifted downward (25 cm) compared with the four others. Just a small change in solute concentration dramatically changed the bulk resistivity of the subsurface, making the ERT measurements very sensitive to the location of the tracer.
View larger version (15K):
[in this window]
[in a new window]
|
FIG. 14. Relative solute concentration profiles of the five best realizations of the combined analysis as they developed with time. The range of all the realizations of solute curves is indicated with gray.
|
|
|
Discussion
|
|---|
The collected geophysical data used in the integrated data fusion approach presented in this study described infiltration of water and solute in a 12-m-deep section of the unsaturated zone. The collected data set had two main advantages compared with the geophysical data used in Binley and Beven (2003): (i) large dynamic changes in the moisture content in the investigated area were expected as a result of a relatively high infiltration rate, i.e., 88.4 mm d–1; and (ii) additional information concerning solute transport was contained in the cross-borehole electrical resistivity data. Despite these promising new features, the presented approach did not prove to be able to identify the subsurface unsaturated parameters unequivocally. To model the observed infiltration patterns satisfactorily, a five-layered model was necessary and this conceptual geologic model resulted in 20 free parameters. The synthetic tests showed that only Ks and n of the top 2 to 3 layers could be constrained, and that somewhat poorly.
We recognize that more realizations and wider ranges might be needed to constrain the parameters better. The proposed sampling approach was very computationally expensive, however, and ways to improve this should also be sought. Additional data or information, such as pressure measurements (Rucker and Ferré, 2004) and correlation of parameters, could readily be incorporated to constrain the data further and possibly minimize the number of realizations needed. Furthermore, the collected data acquisition scheme could be designed more carefully to minimize the amount of data used in the geophysical forward analysis. The geophysical data used in this study were originally collected for sequential data fusion and, as a result, the number of ERT measurements is most likely unnecessarily high. Furthermore, the experiment was designed as a one-dimensional two-step infiltration, where the tracer was added after 4 d. Mass balance considerations revealed that water was being diverted away from the measurement area; however, the consistency of the water loss during the initial 8 d led us to adopt an effective infiltration rate of half the applied irrigation. We feel confident that this assumption is valid, but one could consider estimating the effective infiltration in future analyses to minimize this uncertainty.
Additional assumptions used in the presented analysis that may have had an impact on the obtained results include:
- The definition and internal weighting of the misfit values. The sensitivity of the sought-after hydraulic parameters toward the information contained in the two measurement types (ERT and GPR) may be enhanced by varying the used equations and the weights, given the data types.
- Assumptions regarding the conceptual geologic model, including the one-dimensional flow assumption and the layered model setup, the assumed constant petrophysical parameters, the use of an effective infiltration, a constant dispersivity value for the entire soil column, and the used a priori parameter ranges. Most of these assumptions were made to simplify the conceptual model and minimize the number of free parameters.
- The initial conditions. An initial background moisture content profile was used in the hydrologic simulations. To avoid secondary water fronts resulting from erroneous retention properties, a preinfiltration period could be included in the simulation. This type of simulation, however, introduces additional uncertainties concerning the measured precipitation and assumptions regarding the estimation of the evapotranspiration.
- The use of only 1 d. In order better to describe the dynamic properties of the subsurface, data from multiple days could improve the ill-posedness of the problem. We observed in Fig. 10, however, that data from six additional days did not improve the constraint on the parameters substantially.
- Finally, it is extremely important to note that this analysis was based on the key assumption that the hydrologic model could be expressed using Richards' equation and the van Genuchten parameterization. These models have been developed based on small-scale experiments and may not be valid at the scale on which we were working. Attempts should therefore be made to test additional parameterization approaches.
|
Conclusions
|
|---|
A framework that allows multiple geophysical data sets to be combined using an integrated data fusion approach has been developed. The collected data were used in a forward manner to evaluate a series of plausible hydrologic models. By avoiding the use of geophysical inversion tools, uncertainties connected with inversion, such as artifacts and subjective constraints (e.g., regularization choice) were eliminated. Nonetheless, the proposed methodology still required assumptions regarding choice of petrophysical relationships, conceptual geologic model, investigated parameter ranges, misfit functions, and data error levels. All these variables are connected with uncertainty and represent subjective considerations that may influence the obtained results markedly.
We used electrical resistance and radar travel time data collected using cross-borehole ERT and GPR, respectively, during a forced solute and water infiltration experiment. Results showed that a five-layered model of the subsurface was necessary to capture the most important infiltration patterns through the top 12 m. In this conceptual geologic model, the water and tracer flow was controlled by 20 independent hydraulic parameters. A synthetic evaluation of the proposed methodology for the given setup showed that the data collected only provided information enough to constrain the saturated hydraulic conductivity, Ks, and the van Genuchten parameter n of Layers 1, 2, and to some extent Layer 3, while it was not possible to constrain the residual moisture content, r, or the van Genuchten parameter . A similar result was obtained when real data were applied. No single optimal parameter set appeared to exist that honored the data. The selected parameter ranges of Ks in Layers 1, 2, and 3 were somewhat constrained, however, and the obtained results were found to be consistent with layered changes in grain size distributions collected at a nearby well.
|
ACKNOWLEDGMENTS
|
|---|
We wish to thank the Water Department, Section for Water Quality, at Copenhagen Energy for all their help and for use of their field location at the Arrenæs infiltration plant. Furthermore, Majken Looms thanks the many helpers that assisted her in the field data collection and, in particular, Jannick K. Jensen, who conducted retention experiments on soil samples collected at Arrenæs.
|
REFERENCES
|
|---|
- Alumbaugh, D., P.Y. Chang, L. Paprocki, J.R. Brainard, R.J. Glass, and A. Rautman. 2002. Estimating moisture contents in the vadose zone using cross-borehole ground penetrating radar: A study of accuracy and repeatability. Water Resour. Res. 38(12):1309, doi:10.1029/2001WR000754.
- Archie, G.E. 1942. The electrical resistivity log as an aid in determining some reservoir characteristics. Trans. Am. Inst. Mining Metall. Petrol. Eng. 146:54–62.
- Beven, K., and A. Binley. 1992. The future of distributed models: Model calibration and uncertainty prediction. Hydrol. Processes 6:279–298.[CrossRef]
- Binley, A. 2007. R3: Summary. Available at www.es.lancs.ac.uk/people/amb/Freeware/r3/r3_summary.htm (verified 20 Dec. 2007). Lancaster Univ., Lancaster, UK.
- Binley, A., and K. Beven. 2003. Vadose zone flow model uncertainty as conditioned on geophysical data. Ground Water 41:119–127.[CrossRef][Web of Science][Medline]
- Binley, A., G. Cassiani, R. Middleton, and P. Winship. 2002a. Vadose zone flow model parameterisation using cross-borehole radar and resistivity imaging. J. Hydrol. 267:147–159.[CrossRef]
- Binley, A., P. Winship, R. Middleton, M. Pokar, and J. West. 2001. High-resolution characterisation of vadose zone dynamics using cross-borehole radar. Water Resour. Res. 37:2639–2652.[CrossRef]
- Binley, A., P. Winship, L.J. West, M. Pokar, and R. Middleton. 2002b. Seasonal variation of moisture content in unsaturated sandstone inferred from borehole radar and resistivity profiles. J. Hydrol. 267:160–172.[CrossRef]
- Cassiani, G., and A. Binley. 2005. Modeling unsaturated flow in a layered formation under quasi-steady state conditions using geophysical data constraints. Adv. Water Resour. 28:467–477.[Medline]
- Cassiani, G., R. Deiana, A. Villa, V. Bruno, A. Bagliani, M. Miorali, and N. Fusi. 2006. A water injection experiment in the vadose zone: The use and value of non invasive cross-hole data for model calibration. In P.J. Binning et al. (ed.) Proc. Int. Conf. on Computational Methods in Water Resour., 16th, Copenhagen, Denmark. 19–22 June 2006. Available at proceedings.cmwr-xvi.org (verified 20 Dec. 2007). Tech. Univ. of Denmark, Lyngby.
- Cassiani, G., C. Strobbia, and L. Gallotti. 2004. Vertical radar profiles for the characterization of deep vadose zones. Vadose Zone J. 3:1093–1105.[Abstract/Free Full Text]
- Constable, S.C., R.L. Parker, and C.G. Constable. 1987. Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics 52:289–300.[CrossRef][Web of Science]
- Daily, W., A. Ramirez, D. LaBrecque, and W. Barber. 1995. Electrical resistance tomography experiments at the Oregon Graduate Institute. J. Appl. Geophys. 33:227–237.[CrossRef]
- Daily, W., A. Ramirez, D. LaBrecque, and J. Nitao. 1992. Electrical resistivity tomography of vadose water movement. Water Resour. Res. 28:1429–1442.[CrossRef]
- Day-Lewis, F.D., K. Singha, and A.M. Binley. 2005. Applying petrophysical models to radar travel time and electrical resistivity tomograms: Resolution-dependent limitations. J. Geophys. Res. 110:B08206, doi:10.1029/2004JB003569.
- Ferré, T.P.A., G. von Glinski, and L.A. Ferré. 2003. Monitoring the maximum depth of drainage in response to pumping using borehole ground penetrating radar. Vadose Zone J. 2:511–518.[Abstract/Free Full Text]
- French, H.K., C. Hardbattle, A. Binley, P. Winship, and L. Jakobsen. 2002. Monitoring snowmelt induced unsaturated flow and transport using electrical resistivity tomography. J. Hydrol. 267:273–284.[CrossRef]
- Hansen, T.M., M.C. Looms, and L. Nielsen. 2008. Inferring the subsurface structural covariance model using cross-borehole ground penetrating radar tomography. Vadose Zone J. 8:249–262 (this issue).
- Hubbard, S.S., J.E. Peterson, E.L. Majer, P.T. Zawislanski, K.H. Williams, J. Roberts, and F. Wobber. 1997. Estimation of permeable pathways and water content using tomographic radar data. Leading Edge 16:1623–1628.[Abstract]
- Kowalsky, M.B., S. Finsterle, J. Peterson, S. Hubbard, Y. Rubin, E. Majer, A. Ward, and G. Gee. 2005. Estimation of field-scale soil hydraulic and dielectric parameters through joint inversion of GPR and hydrological data. Water Resour. Res. 41:W11425, doi:10.1029/2005WR004237.
- Kowalsky, M.B., S. Finsterle, and Y. Rubin. 2004. Estimating flow parameter distributions using ground-penetrating radar and hydrological measurements during transient flow in the vadose zone. Adv. Water Resour. 27:583–599.[CrossRef]
- LaBrecque, D.J., A.L. Ramirez, W.D. Daily, A.M. Binley, and S.A. Schima. 1996. ERT monitoring of environmental remediation processes. Meas. Sci. Technol. 7:375–383.[CrossRef][Web of Science]
- Linde, N., A. Binley, A. Tryggvason, L.B. Pedersen, and A. Revil. 2006. Improved hydrogeophysical characterization using joint inversion of crosshole electrical resistance and ground penetrating radar traveltime data. Water Resour. Res. 42:W12404, doi:10.1029/2006WR005131.[CrossRef]
- Looms, M.C., K.H. Jensen, A. Binley, and L. Nielsen. 2008. Monitoring unsaturated flow and transport using cross-borehole geophysical methods. Vadose Zone J. 8:227–237 (this issue).
- McKay, M.D., R.J. Beckman, and W.J. Conover. 1979. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21:239–245.[Medline]
- Moysey, S., R. Knight, and K. Singha. 2006. Relating geophysical and hydrologic properties using field-scale rock physics. In P.J. Binning et al. (ed.) Proc. Int. Conf. on Computational Methods in Water Resour., 16th, Copenhagen, Denmark. 19–22 June 2006. Available at proceedings.cmwr-xvi.org (verified 20 Dec. 2007). Tech. Univ. of Denmark, Lyngby.
- Parkin, G., D. Redman, P. von Bertoldi, and Z. Zhang. 2000. Measurement of soil water content below a wastewater trench using ground-penetrating radar. Water Resour. Res. 36:2147–2154.[CrossRef]
- Reynolds, J.M. 1997. An introduction to applied and environmental geophysics. John Wiley & Sons, New York.
- Rucker, D., and T.P.A. Ferré. 2004. Parameter estimation for soil hydraulic properties using zero-offset borehole radar: Analytical method. Soil Sci. Soc. Am. J. 68:1560–1567.[Abstract/Free Full Text]
- Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. ROSETTA: A computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol. 251:163–176.[CrossRef]
- imnek, J., M. ejna, and M.Th. van Genuchten. 1998. The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat and multiple solutes in variably saturated media. Version 2.0. IGWMC TPS-70. Int. Ground Water Modeling Ctr., Colorado School of Mines, Golden.
- Singha, K., and S.M. Gorelick. 2005. Saline tracer visualized with three-dimensional electrical resistivity tomography: Field-scale spatial moment analysis. Water Resour. Res. 41:W05023, doi:10.1029/2004WR003460.[CrossRef]
- Singha, K., and S.M. Gorelick. 2006. Hydrogeophysical tracking of three-dimensional tracer migration: The concept and application of apparent petrophysical relations. Water Resour. Res. 42:W06422, doi:10.1029/2005WR004568.[CrossRef]
- Singha, K., and S. Moysey. 2006. Accounting for spatially variable resolution in electrical resistivity tomography through field-scale rock–physics relations. Geophysics 71:A25–A28.[CrossRef]
- Topp, G.C., J.L. Davis, and A.P. Annan. 1980. Electromagnetic determination of soil water content: Measurements in coaxial transmission lines. Water Resour. Res. 16:574–582.[CrossRef]
- 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]
- Winship, P., A. Binley, and D. Gomez. 2006. Flow and transport in the unsaturated Sherwood sandstone: Characterization using cross-borehole geophysical methods. Geol. Soc. Spec. Publ. 263:219–231.
- Zelt, C.A., and P.J. Barton. 1998. 3D seismic refraction tomography: A comparison of two methods applied to data from the Faeroe Basin. J. Geophys. Res. 103:7187–7210.[CrossRef]
This article has been cited by other articles:
|
|
|
|
 
S. Lambot, E. Slob, J. Rhebergen, O. Lopera, K. Z. Jadoon, and H. Vereecken
Remote Estimation of the Hydraulic Properties of a Sand Using Full-Waveform Integrated Hydrogeophysical Inversion of Time-Lapse, Off-Ground GPR Data
Vadose Zone J.,
August 11, 2009;
8(3):
743 - 754.
[Abstract]
[Full Text]
[PDF]
|
|
|
|
|
|
|
S. Lambot, A. Binley, E. Slob, and S. Hubbard
Ground Penetrating Radar in Hydrogeophysics
Vadose Zone J.,
February 25, 2008;
7(1):
137 - 139.
[Full Text]
[PDF]
|
|
|
|
|
|
|
K. S. Cordua, M. C. Looms, and L. Nielsen
Accounting for Correlated Data Errors during Inversion of Cross-Borehole Ground Penetrating Radar Data
Vadose Zone J.,
February 25, 2008;
7(1):
263 - 271.
[Abstract]
[Full Text]
[PDF]
|
|
|
TOP
ABSTRACT
INTRODUCTION
