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


     

Published online 13 May 2005
Published in Vadose Zone J 4:300-309 (2005)
DOI: 10.2136/vzj2004.0094
© 2005 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (4)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Ritter, A.
Right arrow Articles by Vanclooster, M.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Ritter, A.
Right arrow Articles by Vanclooster, M.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Ritter, A.
Right arrow Articles by Vanclooster, M.
Related Collections
Right arrow Time Domain Reflectometry, TDR
Right arrow Solute Transport Models

SPECIAL SECTION: ZNS'03 VADOSE ZONE RESEARCH

Using TDR and Inverse Modeling to Characterize Solute Transport in a Layered Agricultural Volcanic Soil

A. Rittera,*, R. Muñoz-Carpenab, C. M. Regaladoa, M. Javauxc and M. Vancloosterc

a Instituto Canario de Investigaciones Agrarias (ICIA), Apdo. 60, 38200, La Laguna, Spain
b Agricultural and Biological Engineering Dep., University of Florida, 101 Frazier Rogers Hall, P.O. Box 110570, Gainesville, FL 32611-0570
c Dep. of Environmental Sciences and Land Use Planning, Unité Génie Rural, Université Catholique de Louvain, Croix du Sud, 2, BP2, B-1348 Louvain-la-Neuve, Belgium
* Corresponding author (aritter{at}icia.es)

Received 21 June 2004.



   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 APPENDIX:
REFERENCES
 
Volcanic soils exhibit particular physical-chemical properties (i.e., strong and stable natural aggregation and high content of variable-charge minerals) that may influence solute transport. To determine if such techniques like TDR and inverse modeling are useful for analyzing solute transport in volcanic soils, we studied the governing transport processes by means of a miscible displacement experiment of Br– in a large undisturbed soil monolith. Bromide resident concentrations at several depths were monitored successfully with TDR technology, while parameters for the convective–dispersive (CDE) and mobile–immobile (MIM) transport models were estimated by inverse modeling. For the relatively high soil moisture conditions, typical of high frequency-irrigation systems that we considered, Br– was found to move slowly by convection–dispersion. Simulations with the CDE and MIM transport models yielded very similar results. Although Br– is generally assumed to behave as a tracer, we found that this anion in our experiment was subject to adsorption at the bottom part of the monolith. This may be explained by the variable-charge nature of the minerals (Fe and Al oxihydroxides) present in this volcanic soil, which exhibited anion exchange when the pH of the soil solution decreased below the zero point of charge.

Abbreviations: BTC, breakthrough curves • CDE, convective–dispersive equation • EC, electrical conductivity • GMCS, global optimization algorithm • MIM, mobile–immobile model • NMS, Nelder–Mead–Simplex • nMSE, normalized mean squared error • TDR, time domain reflectometry


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 APPENDIX:
REFERENCES
 
ALTHOUGH VOLCANIC SOILS occupy only about 1% of the terrestrial surface (FAO, ISRIC, ISSS, 1998), they are very important because they are among the most productive soils of the planet. However, little research has been done on their transport properties (Magesan et al., 2003). In the Canary Islands (Spain), volcanic soils are crucial since they produce 90% of the main export crops, bananas (Musa acuminata Colla) and tomatoes (Lycopersicon esculentum Mill.). These soils exhibit special properties as a result of the strong natural aggregation of particles, the high concentration of Fe and Al oxihydroxides, and the presence of allophanic clays with large specific surface area and water affinity (Moldrup et al., 2003; Regalado et al., 2003). Considering the strong and stable natural aggregation of volcanic soils, the soil liquid phase is often assumed to be divided into two water regions: a mobile (dynamic) phase and an immobile (stagnant) phase associated with the less permeable region of the soil matrix (Mallants et al., 1994). Evidence of the presence of mobile and immobile regions in the soil of this study was reported by Regalado et al. (2003) and Muñoz-Carpena et al. (2005).

The intensive use of agrichemicals in volcanic agricultural fields, as for many other areas of the world, has led to groundwater contamination. Especially in Tenerife (Canary Islands), widespread groundwater pollution has been reported in the traditional agricultural areas along the coast (Muñoz-Carpena et al., 2001). This is a continuing problem that requires the minimization of agrichemical leaching losses. In this context, field-tested numerical leaching models are recommended for studying solute transport through the vadose zone to the underlying groundwater. Such models can be useful tools for understanding the movement of solutes in the soil and for evaluating the potential effect of alternative agricultural practices on groundwater contamination. Performing a representative study of the solute transport requires the application of the simulation model at the field scale, or at the scale of large undisturbed soil cores (monoliths). The use of numerical simulation models further requires the estimation of vadose zone flow and transport parameters, which cannot be measured directly in most cases (Jacques et al., 2002). Also, when considering layered soil profiles, the number of transport parameters may increase considerably.

Miscible displacement experiments are suitable for calibrating solute leaching models since they provide information about such processes as preferential flow, hydrodynamic dispersion, ion exchange, and adsorption under various flow rate and soil water content conditions (Ersahin et al., 2002). These experiments usually involve the application of a solute pulse at the soil surface, followed by measurements of the solute flux and/or resident concentration in the profile. Obtaining solute breakthrough curves (BTCs) from large soil monoliths is not an easy task. Traditional techniques for solute concentration measurements (e.g., soil coring and solution extractors) are usually inappropriate for obtaining high quality data with good spatiotemporal resolution. For this reason, time domain reflectometry (TDR) has become increasingly popular as it allows for continuous and simultaneous measurements of the soil water content ({theta}) and the electrical conductivity (EC) of the soil solution. When the tracer is a saline solute, and for certain temperature conditions and low background salinities, changes in EC can be linearly related to changes in the solute concentration. Time domain reflectometry is a less-destructive and more cost-effective method enabling continuous readings at different soil depths (Vanclooster et al., 1995). The use of TDR for characterizing solute transport has been reported for both laboratory (e.g., Mallants et al., 1994; Heimovaara et al., 1995; Vanclooster et al., 1995; Vanderborght et al., 2000; Seuntjens et al., 2001; Ersahin et al., 2002; Javaux and Vanclooster, 2003) and field studies (e.g., Jacques et al., 1998). However, although it is well known that volcanic soils generally exhibit atypical dielectric responses that affect TDR soil moisture measurements (Tomer et al., 1999; Miyamoto et al., 2001; Regalado et al., 2003), the implications for soil EC determination have received little attention (Vogeler et al., 1996). The soil relative dielectric permittivity is a complex number whose real part accounts for the soil water content and whose imaginary component reflects ionic conductivity losses. Since the EC of the soil solution and the imaginary dielectric constant are interrelated, dielectric peculiarities related to the real part of the permittivity of volcanic soils are expected to also affect the imaginary part (i.e., EC; Muñoz-Carpena et al., 2005). A comprehensive review of advances in dielectric and electrical conductivity measurement in soils using TDR can be found in Robinson et al. (2003).

In addition, there is a need for developing efficient and objective methodologies for estimating solute transport model parameters. An increasingly popular procedure for estimating solute transport parameters is inverse modeling, where an optimization algorithm is coupled with a forward transport model. While inverse modeling of soil hydraulic parameters is common, calibration of solute transport properties using this method is not widespread (Hopmans et al., 2002) because this approach requires a reliable and detailed (in time and/or space) data set of solute transport, which is often difficult to obtain (Jacques et al., 2002).

The classical approach of modeling the transport of a tracer in soils is represented by the CDE (Biggar and Nielsen, 1967), which considers equilibrium transport of a nonsorbing, nonreactive solute in a one-dimensional flow system. In soils with aggregates or large macropores, rapid transport through the mobile phase can cause early solute breakthrough, while diffusion of solute from immobile water back to the mobile region can produce BTC tailing. The resulting BTCs may be asymmetrical and often cannot be described using the CDE. van Genuchten and Wierenga (1976) described an alternative two-domain model known as the mobile–immobile model (MIM). Convective–dispersive transport in this approach is assumed to occur only in the mobile water domain, while water in the immobile region is not available for convective transport, but acts as a source or sink for solutes for the mobile phase. Solute exchange between both regions is diffusion controlled and described by means of a first-order rate exchange process.

In this study we tried to determine whether TDR measurements in combination with inverse modeling are useful to quantify solute transport in a layered volcanic soil with particular physical and chemical properties. Solute transport through a large, layered, agricultural volcanic soil monolith was characterized by conducting a miscible displacement experiment using a Br– (KBr) pulse. Volume-averaged resident concentrations were monitored with TDR at different soil depths. Transport properties were estimated by inverse modeling with a water flow and solute transport numerical model coupled with a global optimization algorithm. The specific objectives of this study were (i) to validate the use of the TDR at relatively high soil moisture conditions to monitor saline solute resident concentrations in volcanic soils at several depths, (ii) to characterize the transport of nonsorbing, nonreactive solutes in this type of soil, and (iii) to evaluate whether the mobile–immobile regions are important for solute transport under relevant field hydraulic conditions (i.e., high water contents).


   MATERIALS AND METHODS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
RESULTS AND DISCUSSION
 CONCLUSIONS
 APPENDIX:
REFERENCES
 
Experimental Set-Up
The solute transport study was conducted in an undisturbed volcanic soil column taken from a banana (‘Giant Cavendish’) field in Tenerife (Canary Islands, Spain). Since the 15th century, terraced fields have been constructed in the coastal areas on top of weathered and fractured basaltic rock (old lava flows) by first building a retaining rock wall on the steep volcano slopes. A bottom drainage layer of crushed basaltic rock (20 cm deep) is subsequently put into place, after which a layer of the transported soil (70–90 cm) is added. Under the traditional production system (surface irrigation), the soil layer was replaced with new soil every 50 yr at the time when the plantation was renewed.

The soil in this study is an Andisol with well-developed andic characteristics (ISSS, ISRIC, FAO, 1994; Moldrup et al., 2003), that is, strong natural microaggregation that translates into high water retention, porosity, specific surface area, and saturated hydraulic conductivity.

To extract the column of undisturbed soil, a custom hydraulic press was used to insert a stainless-steel cylinder (85 cm, 45-cm diam., 0.4-cm wall thickness) slowly into the soil. Once inserted, the cylinder was isolated by excavating the surrounding soil. After covering the top and bottom with appropriate caps, the cylinder was transported to the laboratory. Figure 1 presents a sketch of the laboratory experimental set-up. The soil monolith was equipped with 21 TDR probes (three 20-cm rods of 0.3-cm diam. with a 2.5-cm separation) for measuring the soil water content and the solute concentration at seven depths (denoted as A–G). They were inserted horizontally, 10 cm apart in the vertical direction starting from the top. At each depth, three TDR probes were inserted at 120° from each other. This increased the sampling region, thus ensuring detection of the solute plume, and guaranteed effective one-dimensional concentrations by averaging (Javaux and Vanclooster, 2003). All probes were muliplexed and connected to a TRASE TDR device (Soilmoisture, Inc., Santa Barbara, CA). In addition, two solution extractors (100 mm, 2.5-mm diam.; Rhizon, Eijkelkamp, Giesbeek, the Netherlands) were inserted horizontally into the monolith at each of the seven depths. A suction of 600 cm was applied to the 14 extractors to sample the soil solution periodically. Since the volumes sampled ({approx}20 mL) and the areas of influence of the extractors were small, we believed that they did not significantly affect the fluid flow and the solute resident concentrations. Temperature was monitored with a thermistor inserted into the center of the soil column.



View larger version (27K):
[in this window]
[in a new window]
  		Fig. 1. Experimental set-up for transport experiments in the volcanic soil monolith.

 
The soil monolith was placed on top of a 5-cm-thick saturated sand bed (73-µm particle size), which was connected to a constant-level reservoir using transparent tubing. Thus, by setting the reservoir at some distance vertically from the bottom of the column, while maintaining continuity, a constant suction head could be applied (Fig. 1). Irrigation was applied to the top with a small rainfall simulator that was constructed using a 550- by 550- by 32-mm Plexiglas box equipped with 310 hypodermic needles (6 mm, 0.3-mm diam., spaced 20 mm apart) placed through and glued onto the bottom. The solution was pumped to the rainfall simulator from a large container. Continuous readings of the solution level in the container were used to estimate irrigation flow rates. A collector, equipped with a pressure transducer, was used to continuously measure the volume of water leaving from the base of the monolith (effluent) during the experiment. Custom PC software (developed at I.T.A.C.L., Valladolid, Spain) was used to initiate and log readings (i.e., TDR soil moisture contents, solute concentrations, outflow rates, and temperature) automatically during the experiment.

In a previous study, Ritter et al. (2004) reported that the monolith consisted of four horizons having different water retention properties; they provided soil hydraulic properties of each horizon as obtained by inverse modeling (Table 1).


View this table:
[in this window]
[in a new window]
  		Table 1. van Genuchten–Mualem hydraulic parameters and bulk density for the undisturbed volcanic soil column (Table 1 and p. 135 in Ritter et al., 2004)){dagger}

 
Miscible Displacement Experiment
The miscible displacement experiment was performed in three steps. The monolith was first irrigated with a background solution until the electrical conductivity (measured with a laboratory EC meter) of the soil solution collected in suction samplers at all depths was the same as that of the bottom outflow. A solution of 0.005 M CaSO4 was used to avoid soil dispersion, while thymol was added to serve as a microbial inhibitor (Dane and Hopmans, 2002). In a second step, approximately one pore volume of a 0.025 M KBr tracer solution was applied at a quasi-fixed flow rate of 1.7 ± 0.2 mm h–1 for 250 h. Maintaining the flow rate constant was not possible because of clogging problems in the rainfall simulator. In the third and last step, the irrigation solution was changed again to the background solution for an additional 710 h. The bottom boundary was set at 10-cm suction during the experiment, which is in the range of average field values measured at that depth in a previous study (Muñoz-Carpena, 1999). The initial soil water status is given in Table 2.


View this table:
[in this window]
[in a new window]
  		Table 2. Soil water status during the experiment.

 
Solute flux concentrations were estimated by measuring the EC of samples taken periodically from the effluent. Resident concentrations at the seven observation depths were measured using soil solution samples collected with the suction extractors and also estimated from TDR measurements. The latter approach is based on the assumption that the electrical conductivity of the bulk soil (ECa) and that of the soil solution (ECw) at a particular water content are linearly related for salinity levels between 1 and 50 dS m–1 (Rhoades et al., 1976; Ward et al., 1994). Furthermore, assuming a linear relationship between the electrical conductivity and concentration of the soil solution, the relative concentration at any depth and time can be described by

[1]
where c(z,t) and C(z,t) are the relative and absolute concentrations at depth z and time t, respectively, and subscripts "o" and "i" denote input and initial concentrations, respectively. According to an equation initially proposed by Rhoades et al. (1976), Muñoz-Carpena et al. (2005) found that ECw of this volcanic soil can be estimated (R2 = 0.986) from soil water content, {theta} (L3 L–3) and ECa (dS m–1) measurements as follows:

[2]

Both ECa and {theta} can be measured with TDR. The soil water content was estimated using a specific TDR calibration for the same soil as used in this study (Regalado et al., 2003). According to Nadler et al. (1991), ECa is related to the impedance of electromagnetic wave moving through the soil as follows:

[3]
where Kcc is the cell constant of the TDR probe (m–1), Z is the soil bulk impedance ({Omega}), and ft is a temperature correction factor (ft = 1 at 25°C). To account for cable losses and the presence of connectors, the multiplexer or other discontinuities in the transmission line, we calculated Z as proposed by Castiglione and Shouse (2003):

[4]
where Z0 is the characteristic impedance of the coaxial cable (50 {Omega}), {rho}0, is the sample reflection coefficient, and {rho}air and {rho}sc are the reflection coefficients measured in air and in the short-circuited probe, respectively. Reflection coefficients were calculated as described in Muñoz-Carpena et al. (2005).

The Kcc of each TDR probe was obtained by immersing the probe in six different KBr solutions of known concentration (Heimovaara et al., 1995) ranging from 0.5 to 4.0 dS m–1. The average fitted Kcc was 112 ± 24 m–1. The background electrical conductivity ECw,i at concentration Ci was obtained from the initial conditions before the tracer application (Step 2 of the experiment). On the other hand, since the applied solute pulse was relatively long, the ECw,o at concentration Co was estimated from the TDR readings corresponding to the maximum (saturation) of the BTC (Mallants et al., 1994). At those depths where the input concentration was not reached, ECw,0 was set equal to the EC of the input solution.

The Forward Numerical Model
We selected the mechanistic-deterministic WAVE model (Vanclooster et al., 1996) to describe the flow and transport processes in this soil monolith with different horizons. This code simulates the one-dimensional transport of solute and water in the vadose zone. Transient flow is described with the one-dimensional, isothermal Richards equation for a variably saturated, rigid porous medium, using the mass-conservative numerical scheme proposed by Celia et al. (1990). The soil moisture retention curve is assumed to be of the form given by van Genuchten (1980), while the unsaturated hydraulic conductivity function is described with the van Genuchten–Mualem model (Mualem, 1976; van Genuchten, 1980).

Considering equilibrium (i.e., homogeneity and perfect solute mixing), solute transport of a nonsorbing, nonreactive solute in a one-dimensional flow system reduces to the CDE:

[5]
where C is the solute concentration of the soil solution (M L–3), {theta} is the soil water content (L3 L–3), t is the time (T), z is the vertical distance from the soil surface (L), D is the apparent dispersion coefficient (L2 T–1), and v the average pore water velocity (L T–1) equal to the Darcian flux divided by {theta}. D accounts for both chemical diffusion and hydrodynamic dispersion (Vanclooster et al., 1996). Neglecting chemical diffusion, D can be defined as D = {lambda}v, where {lambda} is the hydrodynamic dispersivity (L).

The MIM for nonsorbing, nonreactive solute transport during transient flow is written in WAVE as follows (Vanclooster et al., 1996):

[6]

[7]
where Eq. [6] describes solute transport in the mobile region and Eq. [7] transport in the immobile domain. The subscripts "m" and "im" indicate the mobile and immobile soil regions, respectively. Dm is the dispersion coefficient in the mobile phase (L2 T–1), and {omega} is the mass-transfer coefficient, which controls the exchange between both regions (T–1). To account for the two liquid phases, the model uses the parameter ß, which expresses the fraction of total water that is mobile, such that {theta}m = ß{theta} and {theta}im = (1 – ß){theta}. Equation [6] and [7] are subject to the following initial and boundary conditions:

[8]
where tp is the duration of the applied solute pulse (T), and L is the solute transport length (L).

Transport Eq. [5] to [8] are solved using a Crank–Nicolson finite difference scheme. The original solution implemented in WAVE was subject to numerical dispersion (Vanderborght et al., 2002, 2004). To reduce this numerical dispersion to negligible values, an empirical correction term was introduced, which was derived by comparing the results of analytical solutions of a series of well-defined transport experiments with the numerical results.

Formulation of the Inverse Optimization Problem
The transport parameters were estimated using inverse modeling by minimizing the following objective function:

[9]
where the right-hand side represents deviations between observed (c*) and predicted (c) space–time concentration using parameter vector b, n is the number of measurements through time at each depth, and nz denotes the number of observation depths (A–G and the monolith's outlet). The weight of a particular measurement, wi, denotes the measurement error and is set equal to {sigma}i–2, where {sigma} is the standard deviation calculated from values obtained with three TDR probes at each observation depth. Model adequacy, uncertainty, and correlation associated with the estimated parameters were determined according to Hollenbeck et al. (2000) and Ritter et al. (2004).

The global optimization algorithm, GMCS, proposed by Huyer and Neumaier (1999) was used here to minimize the objective function. Previous studies (Lambot et al., 2002; Ritter et al., 2003, 2004) showed that the GMCS, combined sequentially with the Nelder–Mead–Simplex (NMS) algorithm (Nelder and Mead, 1965), was a useful tool for estimating the soil hydraulic parameters. The computer code used in this study was initially developed by Lambot et al. (2002) and later modified by Ritter et al. (2004). For the optimization of the solute transport parameters we further modified this code by coupling the GMCS–NMS algorithm with the solute transport module of the WAVE model.

The goodness of fit of the simulations with the optimized parameters was evaluated in terms of the normalized mean squared error (nMSE) (Wilson, 2001), and by visual inspection of the observed and predicted BTCs. The nMSE per range of observed values, which expresses the proportion of the variance about the 1:1 line compared with the variance of the observed data , was calculated as follows:

[10]

Inverse modeling was performed using the CDE model first and then the MIM approach. Since the monolith contained four horizons with different water retention properties (Table 1) as described in Ritter et al. (2004), we assumed different transport parameters for each horizon. For the CDE we optimized the dispersivity, {lambda}, of the four horizons, while for the MIM approach we also included the nonequilibrium parameters ß and {omega}. All other parameters (notably the soil hydraulic parameters listed in Table 1) were kept constant during the inverse analysis.

Breakthrough Curve Data Analysis
An analysis of the shape of the BTC can provide useful information about several variables. These include the mean breakthrough time, {tau} (when the center of mass of the solute front reaches a given depth); the variance, var (spread relative to {tau}); the coefficient of skewness, SK (BTC symmetry); and the average pore water velocity, v. These variables were calculated using the method of temporal moments proposed by Valocchi (1985). However, since irrigation at the top of the monolith was not applied at a constant rate, we analyzed the BTCs in terms of cumulative irrigation moments. This requires that results are expressed in terms of millimeters instead of hours. Thus, we used the average flow rate (1.7 mm h–1) to express results in terms of adjusted time units. A moment analysis of the flux concentrations (as measured in the effluent) was not included because the effluent was sampled too infrequently.


   RESULTS AND DISCUSSION
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
CONCLUSIONS
 APPENDIX:
REFERENCES
 
Soil solution electrical conductivities obtained with the suction extractor samples were found to closely match the predictions from TDR readings. As an example, Fig. 2 shows the BTCs measured with TDR (three probes) and those obtained with the solution extractors at the middle depth in the monolith (Depth D in Fig. 1). Furthermore, Fig. 3 compares the terms on the right- and left-hand sides of Eq. [2] written as (ECa – 0.112)/ECw = 1.876{theta}2 – 0.512{theta}



View larger version (15K):
[in this window]
[in a new window]
  		Fig. 2. Breakthrough curves observed at the middle of the monolith (Depth D) by TDR (lines) and solution extractors (each with a different symbol).

 


View larger version (21K):
[in this window]
[in a new window]
  		Fig. 3. Correlation between values measured in the suction extractor samples and estimations obtained with the TDR.

 
We made this comparison since most ECw values were around ECw,i or around ECw,o, while different combinations of ECa and {theta} may result in the same ECw. Thus, rearranging Eq. [2] instead of plotting predicted vs. measured ECw values is more appropriate to validate the TDR ECw estimates. Using the electrical conductivity data measured with the suction extractor samples (ECw) and the ECa and {theta} values estimated with the TDR produced a satisfactory correlation between the two terms.

These results confirm the applicability of TDR to quickly and nondestructively estimate the electrical conductivity of the soil solution of this volcanic soil. Furthermore, using the approaches of Nadler et al. (1991) and Castiglione and Shouse (2003), in conjunction with the model of Rhoades et al. (1976), proved to be a good method for monitoring Br– transport in our volcanic soil.

Table 2 shows that soil moisture contents estimated with the TDR probes increased from about 0.41 cm3 cm–3 at the top of the monolith to about 0.52 cm3 cm–3 at the bottom. The average values and their standard deviations indicate that the water content remained almost constant throughout the experiment. In addition, the low standard deviations indicate that the horizontal water content distribution was uniform. The RMSEs (see Appendix) for solute resident concentration at the seven observation depths (Table 3) were calculated using the data set measured with the three TDR probes and, as the predictive variable, the average values of the three TDR probes at each depth. The low RMSEs obtained indicate that our analysis based on averaged relative solute concentrations is justified.


View this table:
[in this window]
[in a new window]
  		Table 3. Breakthrough curve RMSEs and results of the moment analysis.{dagger}

 
Figure 4 (symbols) presents average experimental BTCs obtained from TDR readings at the seven observation depths, and the effluent BTC. Table 3 shows the results of the moment analysis. Notice how the difference in breakthrough time ({Delta}{tau}) is variable, being larger at the bottom of the monolith, especially at Depth G. The spreading of the solute pulse around the mean breakthrough time (var) is also greater at the deeper depths. The average pore water velocity (v) was found to decrease exponentially with depth from 7.43 to 0.26 mm h–1. The solute pulse was almost immediately detected at the first depth, but required approximately 7 d to reach the bottom of the monolith (Fig. 4).



View larger version (20K):
[in this window]
[in a new window]
  		Fig. 4. Average breakthrough curves at the seven observation depths and on effluent. Breakthrough curves were obtained from TDR (symbols) and simulated by the WAVE model with CDE (solid lines), and with MIM (dashed lines).

 
The coefficient of skewness (SK) provides a good measure of the shape of a BTC. Relatively symmetrical BTCs (SK values close to zero) were observed at the first six depths (A–F), thus suggesting equilibrium solute transport. Asymmetrical curves generally indicate the presence of some type of nonequilibrium transport process (Mallants et al., 1994). The relatively slow approach to the input concentration at Depth G suggests nonequilibrium transport. Furthermore, Fig. 4 and the {Delta}{tau} values given in Table 3 indicate delayed Br– breakthrough at Depths F and G. This was confirmed later with the inverse procedure when we also optimized the distribution coefficient (kd) at these depths. When a solute is adsorbed and linear, instantaneous and reversible adsorption is considered, kd (L3 M–1) gives the relation between the adsorbed and dissolved concentrations (Vanclooster et al., 1996). Hence, using the CDE approach, {lambda}H1, {lambda}H2, {lambda}H3, {lambda}F, {lambda}G, kdF, and kdG were simultaneously optimized, resulting in dispersivities equal to 45.6 ± 18.1, 32.4 ± 32.7, 41.6 ± 77.1, 62.4 ± 65.7, and 267.2 ± 103.3 mm for horizons H1, H2, H3, and Depths F and G, respectively, while the distribution coefficients at Depths F and G were found to be kdF = 0.61 ± 0.18 cm3 g–1 and kdG = 2.08 ± 0.20 cm3 g–1. Although the optimized dispersivities exhibited large uncertainties, they showed differences along the profile and, except for H4, are within the range for volcanic soils (1–120 mm) as previously reported by several authors (Vogeler et al., 2000; Magesan et al., 2003). In addition, no parameter correlation was found (correlation coefficients < 0.5).

In contrast to our findings, other authors found a constant, decreasing, or even increasing dispersivity with depth according to the flow conditions imposed (Vanclooster et al., 1995; Javaux and Vanclooster, 2003; Magesan et al., 2003). In a study of solute transport in morphologically different Spodosols, Seuntjens et al. (2001) reported differences in the local transport processes among the various soil layers. Our results may be related similarly to the different horizons in the monolith. Because of the relatively large parameter uncertainty the estimated dispersivity values must be considered with some caution. While a comparison with parameters obtained by fitting an analytical solution to the data would be interesting, this is not so straightforward because of the presence of soil layers with nonuniform soil water contents along the profile. Also, the optimization was based on using resident and flux concentration data simultaneously. An attempt was made to optimize dispersivities for each observation depth simultaneously, but similar large uncertainties were observed for the parameters, and no improvement in the forward simulation was achieved. Visual inspection of the simulated and observed BTCs (Fig. 4) and the calculated nMSE (Table 4) shows that the model described the BTCs satisfactorily at all depths when using the optimized parameters. In general, model predictions were better for the first six depths. Solute concentrations at Depth G and of the effluent were described well when the solute front breaks through; however, the elution parts of the curves were predicted poorly. The results furthermore confirm the considerable delay in Br– breakthrough in the lower parts of the monolith, showing especially high retardation factors at Depths F (R = 2.11 ± 0.01) and G (R = 4.78 ± 0.05). Katou et al. (1996), Vogeler et al. (2000), and Magesan et al. (2003) also reported Br– retardation in volcanic soils; however, these studies obtained much lower R values (1.2–1.8).


View this table:
[in this window]
[in a new window]
  		Table 4. nMSE obtained with the forward simulations with the optimized parameters using CDE and MIM.

 
We next used the MIM approach in attempt to optimize simultaneously the dispersivity ({lambda}), the fraction of mobile water (ß), and the mass transfer coefficient ({omega}) at Depth G, where an asymmetrical curve was observed. Values of {lambda} and kd for the other depths were fixed to the values estimated using the CDE approach. This inverse analysis yielded {lambda}G = 141.5 ± 102.6 mm, ßG = 0.282 ± 0.292, and {omega}G = 2.75 x 10–3 ± 2.75 x 10–3 h–1. A comparison between the calculated nMSE for both approaches (Table 4) shows that in general the MIM improved the predictions at Depth F and of the effluent. When considering the entire profile, using MIM instead of CDE produced only a small decrease in nMSE. Figure 4 includes model predictions using the MIM (dashed lines). In general, the MIM results were very similar to those obtained with the CDE. Comegna et al. (2001) also found small differences between the CDE and MIM analysis of chloride BTC from short undisturbed soil columns from different sites in southern Italy (sandy and clayey soils with bulk densities ranging from 1.24 to 1.48 g cm–3). Also, the high dispersivity for the lower part of the soil column fitted as obtained with the CDE may have been an artifact due to possible nonequilibrium and for experimental complications associated with the lower boundary of the monolith.

Our results indicate that for the aggregated soil, the multilayer CDE approach accurately described solute transport for continuous solute applications. Whereas the MIM was expected to perform considerably better, the success of the CDE approach may be explained in terms of the relatively high moisture conditions maintained by continuous irrigation during the experiment. Under these conditions the solute regions may reach an equilibrium state within a short time (i.e., a small characteristic diffusion time into and out of the immobile phase compared with the time required to develop appreciable concentration changes in the mobile phase). This issue was previously analyzed by Vanderborght et al. (1997). They concluded that instantaneous rather than continuous solute application will result in a higher degree of nonequilibrium, and hence more pronounced differences between the CDE and MIM approaches. Some immobile moisture likely was present due to the large amount of water typically being retained in the microaggregates of volcanic soils (e.g., high residual water contents, see Table 1). In summary, for the flow conditions described in this study, most or all of the pore space of the volcanic soil seemed to have contributed to the convective–dispersive transport process. The upper range of soil moisture maintained during the experiment is typically the most relevant in irrigated agricultural and contaminant transport scenarios.

Another interesting issue is the observed delayed arrival of Br– at Depths F and G, and in the effluent. This can be explained by the mineralogy of the volcanic soil, which promotes anion exchange of the clay fraction of the soil. Bromide is usually employed as a tracer and is generally considered to be virtually nonsorbing. However, several studies (Seaman et al., 1995; Katou et al., 1996; Brooks et al., 1998; Vogeler et al., 2000) reported a linear sorption of Br– in soils or sediments that contain considerable amounts of variable-charge minerals. Volcanic soils contain such minerals, usually in the form of Fe and Al oxihydroxides. Regalado et al. (2003) reported a large Fe and Al oxihydroxide content in samples collected during a detailed soil survey of the field where the experimental column was extracted. These minerals have both positive and negative variable charges, depending on the soil solution pH and ionic strength. For pH conditions below the zero point of charge (pHzpc, i.e., the pH where the total charge from cations and anions at the surface is balanced; Seaman et al., 1995), the minerals are positively charged and hence may be subject to anion exchange. The pHzpc for Fe and Al oxihydroxides ranges from 7 to 9 depending on the composition and degree of crystallinity (van Olphen, 1977). Figure 5 shows the pH of the solutions extracted at each depth during the experiment. The pH values at Depths F and G were <7 (pH < pHzpc), which explains the observed retardation in Br– transport.



View larger version (16K):
[in this window]
[in a new window]
  		Fig. 5. pH of the soil solution extracted at the different depths considered. Dashed lines indicate the general range of pHzpc for Fe and Al oxihydroxides (pH = 7–9).

 
Seaman et al. (1995) reported that the number of anion adsorption sites for conditions below the pHzpc will decrease with increasing concentration of the solute. Bromide retardation will hence be reduced when its concentration exceeds the adsorption capacity of the soil matrix. Intensive agricultural practices often are believed to promote the leaching of clay mineral particles through the soil with the irrigation water; these particles then accumulate in the lower part of the profile (Mori et al., 1999). This process would increase the number of anion adsorption sites in deeper parts of the profile and lead to more Br– retardation. This hypothesis should be confirmed based on clay migration studies, which is beyond the scope of this study. Still, our study shows that caution should be taken when using Br– as a tracer for transport in soils with variably charged minerals, especially if the pH of the soil solution is below the zero point of charge of the predominant minerals.


   CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
APPENDIX:
REFERENCES
 
Monitoring of Br– transport through a volcanic soil column using TDR probes at seven observation depths was successful, indicating that the dielectric properties of the volcanic soil did not considerably affect the TDR electrical conductivity estimates. The transport properties of the different horizons were estimated from a miscible displacement experimental data set by inverse modeling using the WAVE model coupled with the GMCS-NMS optimization algorithm. Two approaches for describing the movement of a nonsorbing, nonreactive solute in the soil, the classical CDE and the MIM, were used. For the high soil moisture conditions considered, the CDE described Br– transport satisfactorily, thus suggesting a situation close to convective–dispersive equilibrium transport. Delayed transport of Br– was observed in the bottom part of the monolith. The resulting retardation was related to the high Fe- and Al-oxihydroxides content typical of volcanic soils, including their accumulation at the bottom part of the profile because of illuviation associated with intensive agricultural practices. The presence of these variable-charge minerals near the bottom of the soil column, in conjunction with the low pH of the soil solution (below the zero point of charge), leads to anion-exchange capacity and consequently Br– adsorption. Considering Br– as an inert tracer in such situations may lead to incorrect transport parameter estimation results.


   APPENDIX:
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 APPENDIX:
REFERENCES
 
ROOT MEAN SQUARE ERROR
The root mean squared error (see Lettenmaier and Wood, 1992) was calculated as follows:

[A.1]
where x* and x represent the observed and the predicted variables, respectively; n is the number of measurements through time for each depth; and nk is the number of TDR probes at each depth.


   ACKNOWLEDGMENTS
 
The authors want to thank Dr. J. Álvarez-Benedí (I.T.A.C.L., Valladolid) for developing the data acquisition software and Mr. M. Morawietz (University Friburg) for help in preparing the experimental set-up. This work was financed with funds of the INIA-Plan Nacional de I+D Agrario (Proyecto I+D SC99-024-C2) and Dirección Gral. de Universidades e Investigación de la Consejería de Educación Cultura y Deportes del Gobierno de Canarias. The research was supported by the Florida Agricultural Experiment Station and approved for publication as Journal Series no. R-09686.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 APPENDIX:
REFERENCES
 




This article has been cited by other articles:


Home page
 Soil Sci.Home page
T. Wohling, J. A. Vrugt, and G. F. Barkle
Comparison of Three Multiobjective Optimization Algorithms for Inverse Modeling of Vadose Zone Hydraulic Properties
Soil Sci. Soc. Am. J., January 25, 2008; 72(2): 305 - 319.
[Abstract] [Full Text] [PDF]

Home page
 Vadose Zone JHome page
J. Mertens, R. Stenger, and G. F. Barkle
Multiobjective Inverse Modeling for Soil Parameter Estimation and Model Verification
Vadose Zone J., August 24, 2006; 5(3): 917 - 933.
[Abstract] [Full Text] [PDF]

Home page
 Vadose Zone JHome page
P. F. Germann and D. Hensel
Poiseuille Flow Geometry Inferred from Velocities of Wetting Fronts in Soils
Vadose Zone J., July 26, 2006; 5(3): 867 - 876.
[Abstract] [Full Text] [PDF]

</
This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (4)
Right arrow Citing Articles via Google Scholar
Google Scholar