Published in Vadose Zone Journal 3:1322-1330 (2004)
© 2004 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
ORIGINAL RESEARCH
In Situ Long-Term Chloride Transport through a Layered, Nonsaturated Subsoil. 1. Data Set, Interpolation Methodology, and Results
M. Javauxa,b,* and
M. Vancloostera
a Department of Environmental Sciences and Land use Planning, Université Catholique de Louvain (UCL), Croix du Sud, 2 Bte. 2, B-1348 Louvain-la-Neuve, Belgium
b Currently, Agrosphere Inst., ICG-IV Forschungszentrum Juelich GmbH, D-52425 Juelich, Germany
* Corresponding author (m.javaux{at}fz-juelich.de)
Received 30 September 2003.
 |
ABSTRACT
|
|---|
A long-term data set of Cl– concentrations in an artificial, leaking lake and in the underlying vadose zone was used to assess in situ vertical transport through a tertiary sandy deposit. Since the temporal resolution of the observed Cl– concentrations in the lake was relatively poor, a method was implemented to assess uncertainty on process identification as caused by a poor definition of the boundary condition. For this purpose we performed continuous stochastic simulations of Cl– concentrations in the lake conditioned to the discretely measured Cl– concentrations. This allows one to generate an ensemble of equiprobable time series of top boundary conditions. Each time series was subsequently used in an inverse convective–dispersive (CD) transport model based on the transfer function concept to identify the apparent transport properties of the unsaturated natural porous medium. For each simulation, one optimal parameter set in the least-square sense is then obtained, and the variance of the optimized parameters reflects the uncertainty due to the poor sampling frequency of Cl– in the lake. For the case study presented, we show that the variance component due to poor sampling of the top boundary condition is small for estimating the transport velocities. This component increases with the hydrodynamic dispersivity. However, for this latter parameter, the sampling variance component was still small compared with the total variability between dispersivity at different depths. Physical interpretation of variability and scale dependency of this latter parameter is studied in more detail in an accompanying paper.
Abbreviations: BC, boundary conditions • CD, convective–dispersive • CDE, convection–dispersion equation • cpm, counts per minute • EC, electrical conductivity • PCS, porous cup solution sampler • pdf, probability density function • SC, stochastic–convective • TDR, time domain reflectometry
|
INTRODUCTION
|
|---|
ALTHOUGH CONSIDERABLE INTEREST exists today in studies of solute transport through natural unsaturated porous media, relatively few large-scale in situ transport experiments have been performed that allow one to validate unsaturated solute transport theories at the scale of interest. Indeed, even if the different processes that determine chemical fate and transport in unsaturated porous media are well known, their interactions and the effects of natural heterogeneity on effective large-scale transport are still poorly understood. This is mainly due to our inability to properly characterize medium heterogeneity, which has a significant effect on the flow field and hence on the transport process. Although new promising geophysical measuring techniques such as ground penetrating radar, electrical resistivity tomography, and magnetic resonance electrical impedance tomography should help us in the near future to characterize the natural heterogeneity more precisely, these techniques are currently insufficient to fully characterize the transport processes at the spatial and temporal scales of interest. Therefore, in-situ tracer experiments remain a reference technique to identify and characterize large-scale transport processes.
The first large-scale in situ unsaturated solute transport experiments were performed during the 1970s by Biggar and Nielsen (1976), Wild and Babiker (1976), Van De Pol et al. (1977), and Starr et al. (1978). Several other field-scale transport experiments were performed later (Jury and Flühler, 1992). These tracer experiments consisted of applying a solute tracer at the surface of an agricultural soil and determining its fate and behavior through conventional soil or soil solution sampling. A common problem with these studies is the low sampling frequency in space and time, which does not allow a robust and well-posed identification of transport processes in unsaturated soil. This problem was partially resolved by the development of the time domain reflectometry (TDR) tracing technique, which allows automatic recording and multiplexing of sensors that can detect an ionic tracer, thereby enabling a higher sampling frequency of the ionic tracer in space and time at the field and plot scales (e.g., Jacques et al., 1998). Another alternative consists of applying different types of dyes. Much effort was recently devoted to refining dye tracer analysis to increase the spatial resolution of sampling (Forrer et al., 2000; Flury and Wai, 2003). However, this method still performs poorly in terms of temporal resolution.
In contrast to controlled field studies focusing on transport within a soil root zone, very few studies deal with transport in the subsoil beyond a depth of a classical soil profile. Indeed, accessibility of the soil for installing tracer probes at these depths is limited. In addition, tracer monitoring in deeper soils takes more time and is more difficult to control. Because boundary conditions in these cases are more difficult to control, unsteady-state transport experiments are much easier to implement. However, the monitored tracer data must then be combined with advanced inverse analysis to elucidate the larger-scale transport process. The power of the inverse identification process will thereby depend on the spatiotemporal resolution with which the boundary conditions (i.e., tracer input) and system state variables (i.e., tracer concentrations in soil) have been defined. Interpolation can be used to numerically decrease the time and space lags among data, as well the boundary condition definition as the system state variable in such inverse analyses, but this likely will introduce additional uncertainty into the identification process. The effect of temporal interpolation on the performance of the inversion process and the uncertainty of the transport parameters as inferred from in situ transport experiments have not yet been discussed in the literature.
We present the results of a long-term solute transport study in an unsaturated sandy Tertiary deposit situated below an artificial, leaking lake and above an unconfined aquifer. The governing transport processes at the scale of the formation are identified by means of inverse modeling using the monitored Cl– concentration data collected in the lake and at seven different depths below the lake. To assess the uncertainty in the identification process due to interpolation of the Cl– concentration data in the lake, stochastic simulations were performed. This technique allows us to generate an ensemble of equiprobable time series, which then can be used as a boundary condition for the inverse solute transport model. For each simulation, one optimal parameter set in the least-square sense is obtained, and the scattering of optimized parameters reflects the uncertainty due to the poor sampling frequency of Cl– in the lake.
Our objective is to present the field site, the available experimental database related to the long-term transport experiment, and the identification methodology. Particular attention is devoted to interpolation of monitored time series involving large time steps and the influence of this interpolation on the uncertainty of the identified transport parameters. An accompanying paper (Javaux and Vanclooster, 2004) focuses on interpretation of the transport process itself and relationships between the transport parameters and morphological features within the sandy deposit.
|
THEORY
|
|---|
Inert Solute Transport
Dispersivity and Mixing Regime
The lateral mixing regime is key to understanding the field-scale transport of water and solutes (Flühler et al., 1996). If the soil is conceptualized as an ensemble of parallel flow regions with different velocities in which an ensemble of particles is released, then the lateral mixing regime will determine the time of lateral mixing of particles between the flow domains compared with the mean travel time of the particles. Two asymptotic cases are defined:- If the particles sample different flow domains and hence experience a range of velocities in a relatively short travel time, then the mixing is considered important and the process is convective–dispersive (CD).
- If the flow domain boundaries are impermeable to particles, and hence no transversal mixing occurs, then the process is stochastic–convective (SC).
For three-dimensional solute transport at the field scale, we expect that the mixing regime will change between these two limiting cases as a function of depth and the imposed boundary conditions. The governing mixing regime can be identified by analyzing the mean and the variance of the solute travel time at different depths. Alternatively, the identified one-dimensional solute transport parameters at different times and depths can be used to discriminate the transport process. Thus, even if the convection–dispersion equation (CDE) cannot predict the dispersion process, the two fitted CDE parameters DL and 
L, the longitudinal dispersion coefficient and the velocity, can still be used to discriminate the "mixing regime" within the flow domain (Vanderborght et al., 2001; Javaux and Vanclooster, 2003). Neglecting diffusion, the longitudinal dispersivity 
L (L) is defined as
 | [1] |
If L is constant or decreases with depth, z, then the lateral mixing process is complete. Alternatively, if L increases with z, lateral solute mixing is incomplete (Vanderborght et al., 2001).
Transfer Function Theory
Solute transport can be modeled using transfer function theory (Jury, 1982). With this theory, the solute flux concentration time series at a given soil depth can be calculated by convoluting a soil characteristic solute transfer function ff with an input flux concentration time series Cfin (Jury and Roth, 1990):
 | [2] |
where ff is the travel time probability density function (pdf) or impulse response function, Cf is the solute flux concentration, and 
in is the solute input time. This equation is subject to five important assumptions: (i) the transfer function ff is measurable, (ii) the transport process is linear, (iii) the system is stationary (i.e., ff is constant), (iv) the water flow is steady, and (v) the profile is initially solute free. If input and output concentration time series are available and ff equation is known and dependent on a number of soil-characteristic parameters, Eq. [2] can be used in an inversion procedure. The corresponding equation for predicting resident concentrations is
 | [3] |
with fr* a transfer function corresponding to the time-normalized resident concentration for a narrow input pulse (Vanderborght et al., 1996).
Stochastic Simulations
Consider a series of observations z(ti), i = 1 ... n, modeled as a realization of a stationary random process Z(t) satisfying
 | [4] |
with µ the global mean, and 
(t) a temporally autocorrelated random residual with a zero mean and a temporal correlation expressed by the semivariogram (Matheron, 1963):
 | [5] |
in which is a vector of temporal lags between t and t + . Knowing the semivariogram, an interpolated value between observations at unsampled time t0 can be calculated using the simple kriging algorithm:
 | [6] |
in which wi are the kriging weights chosen such that there is no bias in the prediction and the prediction error variance is minimized. Kriging leads to a single uncertain time series from a set of measurements. The kriging estimate is, in the least-square sense, the best linear and unbiased interpolator. However, for assessing the impact of uncertainty of the interpolation, alternative time series that preserve the distribution and the correlation structure of the reference dataset must be generated. In that case, a stochastic simulation should be performed (Goovaerts, 1997). Conditional simulation incorporates the observations to simulate possible outcomes at unsampled time events. By carrying out a large number of conditional simulations, one can obtain an estimate of the uncertainty related to the interpolation process. The average value at each unsampled time over a very large number of realizations is identical to the simple kriging estimate.
To reduce the uncertainty in the variogram of the reference variable, auxiliary variables can be used, as long as they have a sufficient correlation with the reference variable. For this purpose, the cross covariance function must be defined, in our case for two variables (Wackernagel, 1998):
 | [7a] |
 | [7b] |
Note that Ci,j is different from Cj,i. The cross semivariance can also be defined as (Vauclin et al., 1983)
 | [8] |
To ensure a stationary random process, the seasonality and the trend have to be removed first from the initial time series. This issue will be discussed below in the section Trend and Seasonality Extraction.
|
MATERIALS AND METHODS
|
|---|
Site Description
The experimental site is located in Louvain-la-Neuve, situated 25 km south of Brussels, Belgium. The subsoil of the region consists of a series of irregular Tertiary deposits in discordance with Paleozoic basal complex, covered by a Quaternary loess material. The Tertiary deposits encompass the Landenean (clayey), Bruxellian (sandy), and Tongrean (clayey) formations. Above the basal complex, the lower part of the Bruxellian sand contains an important aquifer. A textural analysis of the same soil type is given in Javaux and Vanclooster (2003).
The relatively newly built city of Louvain-la-Neuve has a completely separate sewage network for rainfall discharge and waste water. The rainfall discharge is collected in an artificial 5.1-ha lake, which acts also as an infiltration basin. This late has been operational since 1986. The bottom of the lake consists in a 10-cm-thick semipermeable loamy deposit covering the unsaturated Tertiary deposit. The vadose zone below the lake is roughly 6 m deep and composed of Bruxellian sand interbedded with three clayey layers. The Bruxellian aquifer defines the bottom of the vadose zone.
To monitor the quality and quantity of the water leaching from the artificial lake to the aquifer, an observation well was installed through the unsaturated zone in the center of the lake. The observation well was equipped with different probes allowing monitoring of the water content, soil water quality, soil water pressure head, and temperature of the vadose zone (Fig. 1)
.
The Experimental Data Set
Lake Data Set
Water quality of the lake was determined monthly as a composite of water samples collected at four different locations within the lake. Major mineralogical elements (Ca2+, Mg2+, K+, Cl–, NO3–, SO4–), dissolved oxygen, electrical conductivity (EC), and pH were determined using standard chemical techniques.
For the purpose of this study, we used the data for Cl–, which is considered to be a conservative tracer. Chloride concentrations were estimated by high performance liquid chromatography and specific ion electrodes. The large step time series of Cl– data was considered as the first source of information Z1 for the cross-covariogram. The dataset consists of 144 observations starting from 19 Mar. 1986, which was also used as the initial time t = 0.
For 3 yr (from 10 July 1996 to 24 Dec. 1999), EC was measured with a 30-min time step using an automatic EC meter. The resulting time series was used as a secondary information source Z2 for the cross-covariogram and contains 31815 observations.
Observation Well Data Set
The observation well was equipped with seven porous cup solution samplers (PCS) located along a vertical transect sampling the unsaturated zone between the lake and the underlying aquifer (Fig. 2)
. The seven depths are 1.00, 1.90, 2.83, 3.89, 4.40, 4.86, and 5.86 m below the bottom of the lake. The porous cups were located 80 cm from the observation well wall, beyond the influence of water moving down the soil–well interface. It is worth noting that only one PCS was available for each depth. Hence, the measurements were not very representative of the complete formation. Soil solution samples were taken monthly by applying a constant suction of 80 kPa (0.8 bars) during 4 d. Major chemical elements of the soil water were analyzed, leading to a long-term, large time step Cl– time series. These time series, which started 17 Jan. 1989 and ended 28 Nov. 2002, contained 760 observations (Fig. 3)
. The significance of the Cl– time series in terms of water and solute transport processes is discussed in the accompanying paper (Javaux and Vanclooster, 2004).
View larger version (33K):
[in this window]
[in a new window]
|
Fig. 2. Soil profile in the unsaturated zone. Dark gray refers to clay layers, light gray to the first loam layer, and white to sand layers. (Left) Water content distribution and standard variation along the year 1990; (right) the position of porous cup samplers (PCS).
|
|
View larger version (42K):
[in this window]
[in a new window]
|
Fig. 3. Chloride concentration time series measured with porous cups at seven depths within the observation well (the x axis is common for all the graphics).
|
|
Next to the observation well, a neutron probe access tube was installed to monitor water contents in the vadose zone. The probe (Nucleotronics, Denmark) was calibrated each year and returned every 3 yr to the manufacturer for recalibration of the electronics. A linear relationship was used to relate the count rate (counts per minute, cpm) to the wet soil density. A second linear calibration curve was used to relate the cpm to water content. Dry soil density was calculated from the wet density by subtracting the moisture volume fraction. Figure 2 depicts the average water content and standard deviation (for 24 measurements) as measured in the unsaturated zone during 1990. Averaged water content values for 1988 and 1989 were assigned in the same standard deviation range.
Water Flow Boundary Conditions
The water level of the lake was monitored with an automated limnimeter. Half-hourly measurements between 10 July 1996 and 24 Dec. 1999 are given in Fig. 4
. The mean lake depth was 4.07 m, with a standard deviation of 9.49 cm.
View larger version (26K):
[in this window]
[in a new window]
|
Fig. 4. Water flow boundary conditions. (Top) Lake water level; (bottom) aquifer water level. Level 0 corresponds to the lake floor.
|
|
The bottom of the observation well was located in the unconfined aquifer. Monthly measured piezometric data were available from 1 June 1977 to 23 Feb. 1999 and represented a time series of 359 observations. Figure 4 shows the water flow boundary conditions (BC). Based on these observations and the water content profiles reported in Fig. 2, we assumed the water flow BC to be at steady state, so that Eq. [2] or [3] can be used. However, this assumption will be further verified in our second paper (Javaux and Vanclooster, 2004).
Profile Layering
A description of the profile is shown in Fig. 2 (right). The profile was inferred from visual observations of the geological logs collected during installation of the survey well. The profile layering was further confirmed by bulk density measurements with a 
-neutron probe in an aluminum access tube close to the survey well. However, because of the neutron probe sampling window and imprecision in the depth of the geological logs, the clay layers boundaries could not been assessed with a precision <5 cm. Notice that the water content profile is influenced mostly by the textural changes. Water accumulated above each clay layer because of the low conductivity of these layers. The increase in water content in the lower part of the profile was due to the aquifer capillary fringe.
Stochastic Simulation of the Solute Top Boundary Condition
Trend and Seasonality Extraction
The time series definition given by Eq. [4] may include a seasonal component S(t), and a mean defined as a deterministic function of time, the trend µ(t). In this case, the time series should be defined as
 | [9] |
Since we prefer to work with a stationary random process, the seasonality and the trend must be removed first from the initial time series. Typically, the trend is a first- or second-order time function, which can be found by classical curve-fitting techniques. Because the time series for the primary variable (Cl– concentration) covers 17 yr, a trend could be defined. However, we were unable to identify a trend for the auxiliary variable (EC) for a short time period.
To remove the remaining seasonal component, we defined this component as the monthly averaged detrended value for each month of the year:
 | [10] |
in which j represents the year, k represents the month, and 
jk represents the trend estimator for year j and month k. Seasonality for the primary and auxiliary variable time series was slightly different (not shown here), which can be explained by the fact that the Z2 time series covers only 2 yr.
Cross-Covariogram Construction
First, an explanatory correlation analysis was performed to study the correlation between Z1 and Z2. Only 17 pairs were sampled exactly at the same time in both time series. Defining the experimental correlation coefficient as
 | [11] |
with Sxy the experimental covariance between Zx and Zy. We found a 
12 value of 0.71 with a high significance of the regression (P < 0.001).
Next, a nested model composed of a Gaussian model (for short time lags) and an exponential model (for larger time lags) was fitted to the experimental cross-covariogram. To construct the latter, time classes of 3-d durations were defined and all data in each class were set to the average class value, which increased the apparent number of pairs considerably and smoothed the cross-covariogram.
The Cl– covariogram was assessed by scaling the cross-covariogram to the initial Cl– data sill. Indeed, one can easily show that if two variables are related by the equation Z2 = a*Z1 + b, the Z2 variance will be equal to the product of a2 by Z1.
Sequential Gaussian Simulations
The conditional simulations were performed using the sequential Gaussian simulation algorithm (Goovaerts, 1997). With this algorithm, a random path is defined in time for visiting each discretized time moment once. At each moment, the simple kriging algorithm is used to determine the parameters of the local Gaussian conditional cumulated distribution function (mean and variance) based on a normal score covariance model. A value is then drawn from this local conditional cumulated distribution function and added to the sample data set. The two last steps are repeated until all time moments are visited. Hence, each of the simulations respects approximately the histogram of the measurement and the covariance function. With this algorithm, conditioning is done on the neighboring observed Cl– concentrations and previously simulated values of Cl– concentration.
Chloride data were simulated for each day. Finally, the simulated data were back-transformed to incorporate the trend and the seasonality into the simulated Cl– concentrations. All the processing was done in Matlab 6.0 (Mathworks Inc., Natick, MA) using the toolbox BMElib (Christakos et al., 2002) for the statistical analysis.
Modeling
Forward Model
Considering the steady-state water flow condition below the lake, transfer function theory is suited to model the transport. We used the CDE model formalism, which describes solute breakthrough in terms of the dispersivity and a mean transport velocity. Since PCS data may be considered either as resident or flux concentration, Eq. [2] or [3] could be used in this analysis. The transfer functions for flux and resident concentrations in a semi-infinite porous medium, subject to a narrow pulse of solute flux concentration through the inlet at z = 0 and t = 0, are, respectively (Jury and Roth, 1990)
 | [12] |
 | [13] |
where t represents the time and z represents the depth. Equations [12] and [13] are defined such that the infinite integrals of these equations are 1. For the initial conditions, we considered that the subsoil did not contain any Cl– before the first measurement. While this assumption is obviously incorrect, it will be shown not have important consequences on the optimizations.
Inversion Technique
Chloride concentration time series extracted from the PCS of the observation well were used for estimating the CDE transport parameters at each observation depth by inverting Eq. [2] and [3], further developed in Eq. [12] and [13]. An inverse solution was obtained for each realization of the stochastic simulation of the lake Cl– concentration time series, that is, for each realization of the BC. Therefore, 1000 inversions were performed for each depth, resulting in a set of 1000 DL and L parameters. We note that, since the medium is heterogeneous and the solute BC is taken as the interpolated Cl– concentration time series within the lake, the two parameters of Eq. [12] and [13] are effective parameters encompassing the heterogeneity of the medium. Therefore, the L obtained by Eq. [1] is also effective.
To ensure that the problem was well posed, we analyzed the objective function in terms of L and L defined as
 | [14] |
where P* and P are the vectors containing, respectively, the observed and simulated Cl– concentration time series, and e is the vector of residuals. As pointed out by Vrugt et al. (2001), the response surface of the objective function reveals the presence of local minima or a well-defined global minimum, and gives information on the parameter sensitivity and parameter correlations. In case the well-posedness was confirmed, automatic inversion was done with a local search Gauss–Newton algorithm.
Since the lake input time series started at t = 0 (19 Mar. 1986) and the PCS Cl– time series started 1035 d later (17 Jan. 1989), we assumed that the initial solute condition did not influence the objective function; the period between Days 0 and 1035 was considered as an initial internal equilibration for the model.
|
RESULTS AND DISCUSSION
|
|---|
Interpolation of Chloride Concentration Time Series
The experimental and modeled cross-covariograms of Z1 and Z2 are given in Fig. 5
. The simulated 1000 time series of Cl– concentrations in the lake and the 144 known data are shown in Fig. 6
. Simulations are, obviously, passing through the known values since they are conditioned to them. The scattering of simulated values is large when known data in the Cl– concentration time series are sparse, as for example in the months preceding time 3000.
View larger version (62K):
[in this window]
[in a new window]
|
Fig. 6. The 1000 co-conditioned stochastic simulations (in blue) representing 1000 equiprobable Cl– concentrations time series. The red crosses stand for the real Cl– concentration time series.
|
|
Uniqueness of the Inverse Problem
To investigate the uniqueness of the inverse solution, we analyzed the shape of the response surface of the objective function (Eq. [14]) for given depths and input time series. An example is given in Fig. 7
. In this case, a well-defined local minimum is present, which guarantees that the simple local search Gauss–Newton algorithm will be sufficient to find the global minimum of the inverse problem.
Tests were done for a number of randomly selected time series and all depths. Each test revealed similar results as shown in Fig. 7, except for the seventh depth. This was due to poor quality of data obtained from the PCS at that depth (see Fig. 3). This time series was therefore omitted in further analyses.
Uncertainty in Parameters
Uncertainty
At each depth, the 1000 input concentration time series generated by stochastic simulations were used to estimate 1000 optimal values of L and L. Histograms of those 1000 parameters allow visualization of the uncertainty associated with these parameters (Fig. 8)
. Notice that the L distribution is rather narrow. Since the L is related to the first moment of the transfer function, this means that this moment was not altered by the interpolation procedure. Hence, the information contained in the initial time series was sufficient to describe the solute travel time and, thus, L. At each depth, the parameter corresponding to the smallest error between simulation and measurement and the parameter obtained from the simple kriged time series are also reported. One may observe that parameters obtained from simple kriging are located mostly near the peak of the parameter distribution, whereas the optimized velocities with the smallest errors are generally larger. However, given the uncertainty of the parameters derived from a single time series, and due to errors between measured and modeled concentrations, the two optimized velocities are statistically indistinguishable. Moreover, the uncertainty on parameter estimation for one input time series is comparable to the spreading of the parameters for several equiprobable input times. One may conclude that, in this case, the uncertainty of the input concentration is reflected in the error of the model fit. This conclusion, however, cannot be generalized to other cases.
View larger version (38K):
[in this window]
[in a new window]
|
Fig. 8. Histogram of optimized longitudinal velocity values for the seven depths, considering that porous cup sampler concentrations are resident. Red triangles represent parameters obtained with the least error among the 1000 simulations; green triangles represent parameters obtained with the simple kriged input time series.
|
|
Figure 9
shows the histograms of the optimized L. In this case, the uncertainty related to the interpolation is stronger. Larger dispersivity values seem to result in larger uncertainty. However, this uncertainty does not overshadow the strong depth dependency of the dispersivity (Javaux and Vanclooster, 2004). Results indicate that the optimized dispersivities resulting from the kriged time series are generally lower than the averaged dispersivity of the 1000 input time series. Since kriging represents a mean interpolation process, the upper time series will be smoother than for the simulations. Therefore, lower dispersivities are needed to simulated PCS Cl– time series.
View larger version (45K):
[in this window]
[in a new window]
|
Fig. 9. Histogram of optimized longitudinal dispersivity values for the seven depths, considering that porous cup sampler concentrations are resident. Red triangles represent parameters obtained with the least error among the 1000 simulations; green triangles represent parameters obtained with the simple kriged input time series.
|
|
Actually, it is important to note that the dispersivity may be influenced by the interpolation. Indeed, the decrease in the semivariance of the interpolation of the input time series process, or the increase in L of the transfer function, may lead to similar changes in predicted transport, that is, a smoothing of the concentration time series. Both the interpolation procedure and the transfer parameters will affect the variance of the simulated output concentration. This also justifies carrying out an uncertainty analysis on the interpolation if solute transport parameters need to be assessed.
Correlation between L and L
Correlation diagrams characterize the experimental correlation between the two optimized parameters. This correlation may be due to the simulation procedure, which could artificially create a link between the solute arrival times (and therefore the L) and the dispersion (and thus L). It can also be due to the inverse problem itself, although no significant correlation was visible in the calculated response surfaces (Fig. 7). Figure 10
shows the correlation diagrams for each depth, except for the last one. The structure of the correlation is different for each depth, although we could separate them into three groups by comparing their correlation coefficient (): no correlation, < 0.2 (as at the 1-m and 3.89-m depths); a slight positive correlation, 0.2 < < 0.5 (at the 2.83- and 4.4-m depths); and more pronounced positive correlation, < 0.5 at the 1- and 3.89-m depths. Since the correlation structures are depth dependent, it can be presumed that this correlation is not an artifact of the simulation procedure.
Minimum Errors
Each parameter set for different depths and different realizations of the top BC was obtained by inverse modeling. Hence, each parameter set can also be characterized by the residual between measured and simulated concentration in the vadose zone. For each depth, we identified the BC time series which resulted in this smallest residual concentration. In Fig. 8 and 9, triangles indicate the parameter set that produced the smallest residual concentration. There is no physical reason for how these "best" values should correspond to the mean of the calculated distributions, yet it can be noted that the depth dependence of these best parameters is preserved.
Surprisingly, the same simulated input time series resulted in the best (in the least-square sense) travel time series simulated at different depths. This could mean that, even if they are all equiprobable, some simulations are more "realistic" than others. As can be seen in Fig. 11
, the best input Cl– concentration simulations are quite similar and located rather close to the mean simulated value (which is, by definition, equal to the simple interpolation). This confirms the fact that the uncertainty derived from this method is not underestimated. Moreover, it was observed that the same input time series resulted in the minimum error when considering the concentrations as either flux type or resident type.
View larger version (63K):
[in this window]
[in a new window]
|
Fig. 11. Optimal interpolated input time series in green compared with the 1000 input time series generated by stochastic simulation.
|
|
|
CONCLUSION
|
|---|
From 1987 onwards, Cl– concentrations were monitored in an artificial, leaking lake and in the underlying unsaturated sedimentary sandy deposit up to a depth of 5.9 m. This unique long-term data set was used to assess the prevailing solute transport processes in the vadose zone between the lake and the underlying aquifer. Since the temporal resolution of the experimental data set was low, specific attention was paid to interpolation of the solute input time series, which needs to be defined as a continuous input function in forward solute transport simulation. Stochastic simulations conditioned to the measured lake Cl– concentrations were performed and used to assess the parameter uncertainty due to this interpolation. This technique produced an ensemble of equiprobable time series that was used as the upper boundary condition for an inverse solute transport model. Each simulation yielded an optimal parameter set (in the least-square sense), while the variance of the optimized parameters reflected the uncertainty on these parameters due to the undersampled time series used as the top boundary condition.
For our example, this sampling variance component was small for estimation of the transport velocity, but more important for the dispersivity. However, for the latter, the sampling variance component was still small compared with the total observed variability of L. Also, the parameter scattering resulting from the 1000 simulations was comparable to the uncertainty due to the error between simulated and measured soil Cl– time series. However, this cannot be generalized to other cases. It is also important to note that the dispersivity itself may be influenced by the interpolation method. Therefore, we recommend testing the validity of the inverse modeling results by performing an uncertainty analysis each time an interpolation procedure is performed. Variability and scale dependency of the dispersivity will be studied in more detail in the accompanying paper (Javaux and Vanclooster, 2004).
|
REFERENCES
|
|---|
- Biggar, J.W., and D.R. Nielsen. 1976. Spatial variability of the leaching characteristics of a field soil. Water Resour. Res. 12:78–84.
- Christakos, G., P. Bogaert, and M. Serre. 2002. Temporal GIS: Advanced functions for field-based applications. Springer-Verlag, Berlin.
- Flühler, H., W. Durner, and M. Flury. 1996. Lateral solute mixing processes—A key for understanding field-scale transport of water and solutes. Geoderma 70:165–183.[Web of Science]
- Flury, M., and N.N. Wai. 2003. Dyes as tracers for vadose zone hydrology. Rev. Geophys. 41(1):1002. doi:10.1029/2001RG000109
- Forrer, I., A. Papritz, R. Kasteel, H. Flühler, and D. Luca. 2000. Quantifying dye tracers in soil profiles by image processing. Eur. J. Soil Sci. 51:313–322.
- Goovaerts, P. 1997. Geostatistics for natural resources evaluation. Oxford University Press, New York.
- Jacques, D., D.J. Kim, J. Diels, J. Vanderborght, H. Vereecken, and J. Feyen. 1998. Analysis of steady state chloride transport through two heterogeneous field soils. Water Resour. Res. 34:2539–2550.
- Javaux, M., and M. Vanclooster. 2003. Scale- and rate-dependent solute transport within an unsaturated sandy monolith. Soil Sci. Soc. Am. J. 67:1334–1343.[Abstract/Free Full Text]
- Javaux, M., and M. Vanclooster. 2004. In situ long-term chloride transport through a layered, nonsaturated subsoil. 2. Effect of layering on solute transport processes. Available at www.vadosezonejournal.org. Vadose Zone J. 3:1331–1339 (this issue).[Abstract/Free Full Text]
- Jury, W.A. 1982. Simulation of solute transport using a transfer function model. Water Resour. Res. 18:363–368.
- Jury, W.A., and H. Flühler. 1992. Transport of chemicals through soil—Mechanisms, models, and field applications. Adv. Agron. 47:141–201.
- Jury, W.A., and K. Roth. 1990. Transfer functions and solute movement through soil: Theory and applications. Birkhäuser, Basel, Switzerland.
- Matheron, G. 1963. Principles of geostatistics. Econ. Geol 58:1246–1266.[Abstract/Free Full Text]
- Starr, J.L., H.C. DeRoo, C.R. Frink, and J.P. Parlange. 1978. Leaching characteristics of a layered field soil. Soil Sci. Soc. Am. J. 42:386–391.[Abstract/Free Full Text]
- Van De Pol, R.M., P.J. Wierenga, and D.R. Nielsen. 1977. Solute movement in a field soil. Soil Sci. Soc. Am. J. 41:10–13.[Abstract/Free Full Text]
- Vanderborght, J., M. Vanclooster, D. Mallants, J. Diels, and J. Feyen. 1996. Determining convective lognormal solute transport parameters from resident concentration data. Soil Sci. Soc. Am. J. 60:1306–1317.[Abstract/Free Full Text]
- Vanderborght, J., M. Vanclooster, A. Timmerman, P. Seuntjens, D. Mallants, D.J. Kim, D. Jacques, L. Hubrechts, C. Gonzalez, J. Feyen, J. Diels, and J. Deckers. 2001. Overview of inert tracer experiment in key Belgian soil types: Relation between transport, and soil morphological and hydraulic properties. Water Resour. Res. 37:2873–2888.
- Vauclin, M., S.R. Vieira, G. Vachaud, and D.R. Nielsen. 1983. The use of cokriging with limited field soil observations. Soil Sci. Soc. Am. J. 47:175–184.[Abstract/Free Full Text]
- Vrugt, J.A., W. Bouten, and A.H. Weerts. 2001. Information content of data for identifying soil hydraulic parameters from outflow experiments. Soil Sci. Soc. Am. J. 65:19–27.[Abstract/Free Full Text]
- Wackernagel, H. 1998. Multivariate geostatistics: An introduction with application. Springer-Verlag, Germany.
- Wild, A., and L.A. Babiker. 1976. The asymmetric leaching pattern of nitrate and chloride in a loamy sand under field conditions. J. Soil Sci. 27:460–466.
This article has been cited by other articles:
|
|
|
|
 
J. Vanderborght and H. Vereecken
Review of Dispersivities for Transport Modeling in Soils
Vadose Zone J.,
January 24, 2007;
6(1):
29 - 52.
[Abstract]
[Full Text]
[PDF]
|
|
|
|
|
|
|
M. Javaux, J. Vanderborght, R. Kasteel, and M. Vanclooster
Three-Dimensional Modeling of the Scale- and Flow Rate-Dependency of Dispersion in a Heterogeneous Unsaturated Sandy Monolith
Vadose Zone J.,
April 27, 2006;
5(2):
515 - 528.
[Abstract]
[Full Text]
[PDF]
|
|
|
|
|
|
|
M. F. Dyck, R. G. Kachanoski, and E. de Jong
Spatial Variability of Long-Term Chloride Transport under Semiarid Conditions: Pedon Scale
Vadose Zone J.,
September 12, 2005;
4(4):
915 - 923.
[Abstract]
[Full Text]
[PDF]
|
|
|
|
|
|
|
M. Javaux and M. Vanclooster
In Situ Long-Term Chloride Transport through a Layered, Nonsaturated Subsoil. 2. Effect of Layering on Solute Transport Processes
Vadose Zone J.,
November 1, 2004;
3(4):
1331 - 1339.
[Abstract]
[Full Text]
[PDF]
|
|
|
TOP
ABSTRACT
INTRODUCTION
