Published online 27 April 2006
Published in Vadose Zone J 5:529-538 (2006)
DOI: 10.2136/vzj2005.0085
© 2006 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
ORIGINAL RESEARCH
Interpretation of Dye Transport in a Macroscopically Heterogeneous, Unsaturated Subsoil with a One-Dimensional Model
M. Javauxa,b,*,
R. Kasteelb,
J. Vanderborghtb and
M. Vancloostera
a Dep. 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 Agrosphere, ICG-IV, Forschungszentrum Jülich, D-52425 Jülich, Germany
* Corresponding author (m.javaux{at}fz-juelich.de)
Received 12 July 2005.
 |
ABSTRACT
|
|---|
The objective of this study was to identify physical and chemical processes affecting Brilliant Blue (BB) transport in an unsaturated, macroscopically heterogeneous subsoil. We performed a BB leaching experiment in a 1-m-long undisturbed sandy monolith with a 10-cm-thick discontinuous clay layer 0.2 m below the surface. Two-dimensional BB concentration distributions were derived from image analysis at several depths in the monolith. Results showed several features of nonideal transport: (i) a deeper than expected travel depth, (ii) extensive tailing and a secondary peak in the BB depth profile, (iii) a lower than expected BB concentration in the upper part of the soil profile where the invading dye tracer solution was assumed to have replaced the initial solution, and (iv) a decrease in the concentration variability between the surface and the 2-cm depth. These results were in sharp contrast with inert transport experiments that suggested homogeneous flow and transport characterized by small dispersivity lengths. The macrostructure (i.e., the discontinuous clay layer and the nonuniform irrigation) determined the BB transport. A simulation with a one-dimensional transport model using an effective water content and sorption isotherm that was directly derived from the cross-sectional area of the clay layer suggested that the macrostructure alone was not sufficient to explain the BB distribution. Moreover, a two-site kinetic sorption model was needed to reproduce the tailing of the BB profile and the lower than expected concentrations at the soil surface. Since batch sorption experiments did not reveal rate-limited sorption, we postulate that the unsaturated flow conditions led to reduced accessibility of the sorption sites, thus decreasing the sorption rate.
Abbreviations: BB, Brilliant Blue • BTC, breakthrough curve • CDE, convection–dispersion equation • EC, electrical conductivity • TDR, time domain reflectometry
|
INTRODUCTION
|
|---|
BRILLIANT BLUE and dye tracers in general are widely used to visualize flow and transport paths in unsaturated soils. Flury et al. (1994), for example, performed BB leaching experiments in a series of undisturbed soils and showed that structured soils were more prone to bypassing flows than unstructured soils. Numerous other studies have used BB or other tracers to assess preferential flow patterns in soils (see the review of Flury and Wai, 2003). The success of dye tracing in soil hydrology is primarily due to its low costs, the fact that it can be applied on a large scale, and that it gives an almost direct image of the dye distribution with high resolution.
However, dye tracers sorb to the soil matrix, which is a potential drawback when the tracers are used to study water flow. In addition, dye tracers are mostly applied at high concentrations, which further complicates the interpretation of tracer displacement because their sorption isotherms are nonlinear (i.e., concentration dependent). On the other hand, dye tracers may be ideal tools to gain fundamental insights into the transport of nonlinearly sorbing substances (e.g., pesticides) in naturally structured soils.
Few studies have focused on numerical simulations of observed BB patterns. Some studies attempted to model binarized (stained vs. unstained pixels) dye images. Schwartz et al. (1999) reproduced the observed vertical distribution of the stained soil fraction with a random walk model parameterized from batch sorption experiments and inert tracer experiments. They could only describe their observations after increasing the dispersion coefficients, which were derived from the inert tracer experiment, by one order of magnitude. Persson et al. (2001) could satisfactorily model the observed mean and variance of the dye penetration depth in one plot with a diffusion limited aggregation model, of which the parameters (the walking probabilities Pu and Pd) were first fitted to observations from another plot. Olsson et al. (2002) fitted a "universal multifractal" model (Schertzer and Lovejoy, 1987) to the spatial distribution of maximum dye infiltration depths obtained from digitized high-resolution dye images.
Forrer et al. (2000) presented a new procedure to derive high-resolution dye concentration images from dye tracer images. Kasteel et al. (2002) made an attempt to properly reproduce the dye tracer and bromide profiles in a layered soil with a one-dimensional model. They suggested that BB and bromide did not necessarily follow the same flow pathways and thus that BB was not suitable to trace the travel path of water itself; however, as qualitative information, high resolution dye tracer patterns may help to characterize nonlinear transport processes in heterogeneous flow fields, especially when combined with results from other transport experiments.
Numerous modeling studies have considered nonlinear sorbing contaminant transport in heterogeneous porous media (e.g., Liedl and Ptak, 2003; Mansell et al., 1995; Seuntjens, 2002; Young and Ball, 1999). For practical applications (e.g., pesticide registration), simpler models that assume a homogeneous flow field are often used to predict the transport of nonlinearly sorbing components in real soil (e.g., FOCUS Groundwater Scenarios Workgroup, 2000). These models are frequently parameterized based on inert tracer experiments to determine hydrodynamic transport parameters and batch experiments to determine sorption isotherms and kinetics. Nevertheless, this approach may fail for several reasons. First, for the case of nonlinearly sorbing substances, spatially variable local concentrations lead to an additional variability in local retardation constants. Therefore, the spatial averaging of local concentrations in a horizontal plane used in one-dimensional models (homogenization in mathematical terms) may be inappropriate for predicting nonlinear transport in heterogeneous flow fields (Vanderborght et al., 2006). Second, for sorbing substances, the effect of spatially variable sorption constants may have an important impact on transport of reactive substances, which obviously cannot be derived from an inert tracer experiment (e.g., Johnson et al., 2003; Vanderborght et al., 2002) Third, the conditions in batch experiments may differ considerably from conditions in unsaturated and naturally structured soil (e.g., Albrecht et al., 2003; Haws et al., 2004; Sonon and Schwab, 2004). The use of homogenized sieved and dried soil material (homogenization in physical terms) employed for physicochemical sediment characterization may be questioned. Also the impact of nonsaturated conditions on sorption has been investigated by several authors, generally focusing on changes in the retardation factors. Some studies showed positive correlation between retardation and water content (Gamerdinger et al., 2001), while others found the opposite (Estrella et al., 1993; Maraqa et al., 1999).
Besides investigating flow field heterogeneity, which so far has been the main objective of dye tracer studies in soils, distributions of a dye tracer concentration may potentially be used to investigate the effect of chemical heterogeneity, of nonlinear sorption in heterogeneous flow fields, and of sorption in an undisturbed and unsaturated soil matrix on transport.
The objective of our study was to illustrate the potential of dye tracer studies for investigating relevant physical and chemical processes affecting the transport of a reactive substance in an undisturbed soil monolith. For this purpose a BB leaching experiment and a series of inert tracer experiments were performed in an unsaturated heterogeneous undisturbed monolith (0.5 m3). Results of the inert tracer experiments were discussed in Javaux and Vanclooster (2003). The interpretation of BB concentration patterns in horizontal cross sections and the identification of relevant processes are based on comparisons between experimental data and simulations using differently parameterized one-dimensional transport models.
|
MATERIALS AND METHODS
|
|---|
Subsoil
The experiment was performed in a cylindrical undisturbed subsoil monolith that was contained in a large PVC tube (1.05-m height, 0.77-m i.d.). The monolith was extracted 10 m below the ground surface in a sand quarry in Mont-Saint-Guibert (Belgium) using methods described by Vanclooster et al. (1995). The subsoil of the region is a Tertiary marine sandy deposit called the Bruxellean sand, which varies in thickness between 0 and 60 m. Since this deposit was formed by a sedimentation process, it exhibits a series of small horizontal layers (<5 cm thick). Glauconitic clay and sand concretions together with trace fossils are locally present in the sand.
In the laboratory, the soil column was put on top of a 5-cm reservoir filled with fine gravel and with a central outlet gate that collected the outflow. The monolith was equipped with 12 time domain reflectometry (TDR) probes horizontally inserted in three vertical transects (three TDR probes at the 0.15-, 0.45-, and 0.75-m depths, and one TDR probe at the 0.30-, 0.60-, and 0.9-m depths), four horizontally inserted tensiometers (at the 0.10-, 0.30-, 0.60-, and 0.95-m depths), four temperature probes (at the 0.15-, 0.45-, 0.75-, and 0.90-m depths), and one electrical conductivity (EC) meter and one flow meter at the outlet. A detailed description of the equipment is given in Javaux and Vanclooster (2003, their Fig. 1
).
Profile Description
Grain-size distributions were derived from 100-cm3 disturbed samples taken in the monolith sampling trench. The weight fractures of the texture classes ranged from 93.9 to 96.4% for sand, 0 to 1.7% for silt, and 2 to 5.2% for clay (Javaux, 2004).
However, within this rather homogeneous sandy body, macroscopic heterogeneities were observed during the slicing of the monolith after the leaching experiments. The main structural feature was a discontinuous clay layer between approximately 0.20 to 0.30 m below the top surface. The layer was composed of 60% clay and 40% silt. Mineralogic analysis of the clay fraction indicated the presence of smectite and glauconite (or illite), with a small amount of kaolinite. Some sandy concretions and fossils were also found in the lower part of the monolith (around the 0.73-m depth).
Previous Chloride Breakthrough Curves Experiments
Breakthrough experiments were previously conducted using an inert tracer, calcium chloride, under unsaturated and saturated steady-state flow conditions for 11 different flow rates Jw, ranging from 0.02 to 1 times the saturated hydraulic conductivity. Solute breakthrough curves (BTCs) were monitored in the monolith with TDR and in the outflow with the EC meter. Flow and depth dependent apparent velocity and dispersivity parameters were obtained by fitting the solution of the convection–dispersion equation (CDE) to the observed BTCs. A detailed description of these experiments and results is given by Javaux and Vanclooster (2003).
Dye Tracer Experiment
One dye tracer infiltration experiment was performed in the monolith. The organic molecule Brilliant Blue FCF (Food Blue 2, C.I. 42090) was used as dye (for information about the chemistry and adsorption of BB, see Flury and Flühler [1995] and Ketelsen and Meyer-Windel [1999]). First, a constant flow rate of 500 µS cm–1 CaCl2 solution was applied to reach a steady-state water flow of 22.75 cm d–1 corresponding to 50% of the saturated hydraulic conductivity of the monolith. Subsequently, a 1.047 g L–1 BB solution was applied with the same flow rate during 19 h, corresponding to an irrigation dose of 180 mm. After irrigation ceased, the lysimeter was drained by gravity for 8 h, with a zero tension bottom boundary condition.
To control the incoming flow rate during the experiments, a specific irrigation system fed by an accurate pulse pump was designed. It was made of a 0.8 by 0.8 m square reservoir of 1-cm height, with 280 needles (i.d. 0.5-mm) in a 4 by 4 cm grid at the bottom surface of the reservoir. The coefficient of uniformity of the irrigation was 87.5%. The irrigation plate was placed on a 30-cm-high PVC cylinder that was put on the soil monolith and sealed to prevent evaporation during the experiments. The soil surface was covered by a thin layer of gravel to enhance the horizontal redistribution of water dropping from the needles before it infiltrated in the soil.
The monolith was then cut into 17 horizontal slices (at depths of 0, 2, 7.5, 12.5, 18, 20.7, 23, 28.5, 34, 39, 44.4, 50, 55.3, 60.8, 68.4, 73.3, 80.7, and 90.3 cm), and top view pictures were taken at each depth. The soil surface was artificially illuminated by two opposite spotlights of 2 kW focusing at a 45° angle to the soil surface. A wood frame with a ruler was fixed around the circular PVC monolith holder to allow picture corrections for geometrical distortion. Kodack gray and color scales were attached to the panel for color corrections. The slices were photographed with a tungsten color slide film (Fujifilm, Fujichrome 64T, type II, tungsten). At each depth, a picture of a reference green board was taken to correct for inhomogeneous illumination. Kopecky ring soil samples (100 cm3) were taken at 0.15-m depth intervals for determination of soil bulk density and water contents, and additional pictures were taken of Kopecky sampling locations (Javaux and Vanclooster, 2006).
Brilliant Blue Extraction
To relate color and concentration, a few grams of homogeneously stained soil were taken at several locations across the soil surface at each depth. In total, 118 samples were taken. The sample locations were marked with pins and photographed. They were carefully chosen to cover a color range as wide as possible at each depth. In the laboratory, 0.03 to 0.57 g of air-dried soil was filled into a disposable extraction column (8 mL Bakerbond, fritted glass). The BB was extracted from each sample using the same extraction setup as described in Forrer (1997). The samples were consecutively extracted with 2 mL of a 4:1 (% [v/v]) water–acetone solution four to eight times (dependent on initial weight of the soil sample and concentration). Extracts in which BB was visible were either filled up to 20 or to 50 mL (depending on the BB concentration). Extracts without any visible BB were not further diluted. Before measurement, the extracts were filtered through a hydrophilic filter with pore size 0.45 µm. The BB concentration in the extract was derived from UV absorption measurements at 630 nm wavelength using an UVIKON UV spectrometer (calibration in the range of 0.1 to 12 mg L–1) and was afterward recalculated to a BB concentration per unit mass of dry soil.
Image Processing
The procedure proposed by Forrer et al. (2000) was used to process the images. The pictures were first digitized and then corrected for geometrical distortion using known tie points from the ruler. Inhomogeneous illumination was corrected by a background subtraction method. The colors were corrected by adjusting the R, G, and B values to common norms based on the pinned Kodack gray scale. The final pixel size was about 1 by 1 mm. All the image processing steps were performed with Matlab programs.
The calibration data set consisted of 118 samples with BB log-transformed concentration (total concentration per mass of dry soil) and corresponding explanatory variables: red (
), green (
), blue (
), hue (
), saturation (
), value (
), and depth (
). The depth was included as an explanatory variable since the average illumination and background color could have changed with depth. Although HSV variables are somewhat redundant with RGB variables, they were added to the variable set because they could lead to better prediction. We first standardized the explanatory variables to zero mean and unit variance. Subsequently, the Matlab procedure stepwise was used to fit a second-order polynomial of those seven explanatory variables to the data and nonsignificant terms (i.e., terms nonsignificantly different from 0 with a significance level of 95%) were removed. The optimal model was then applied to each pixel to derive concentration maps. The BB concentrations, which were obtained using the regression model and represent mass BB per unit mass of dry soil, were recalculated to BB concentrations defined as mass BB per unit volume of bulk soil using the soil bulk density, which was assumed to be constant in the soil monolith. To check the regression model and the calculated BB concentrations, a mass balance was performed. The total BB mass in the monolith was obtained by integrating the estimated BB concentrations images over depth.
Adsorption Isotherms
Batch studies were performed with BB concentrations of 1, 10, 50, 100, 500, 1000, 2000, and 4000 mg L–1. The dye was dissolved in deionized water as done for the dye tracer experiment in the soil monolith. Duplicate samples of 25 g of air-dried soil and 25 mL of BB solution were used. The soil samples and solutions were mixed in 250-mL plastic centrifugation flasks, shaken for 24 h to equilibrate, and then centrifuged at 14 000 g for 15 min. An aliquot of the supernatant was filtrated through a 0.45-µm filter. The solution absorbency was measured at the 630-nm wavelength in a UV photospectrometer. The relationship between absorbance and BB concentration was calibrated with a series of known BB concentrations.
We could adequately describe the sorption isotherm using the Langmuir equation (Hinz, 2001):
 | [1] |
where 
Cwr is the BB concentration in the liquid phase (mg L–1), S is the sorbed concentration (mass of dye adsorbed per mass of soil) (mgBB kgsoil–1), S
is the maximum or asymptotic sorbed concentration, and is the overall affinity coefficient (L mg–1). Equation [1] was fitted to the loge–transformed data and resulted in the following parameters: S = 1.93 x 103 mgBB kg–1soil and = 1.77 x 10–3 L mg–1. These parameters corresponded well with the predictions of Ketelsen and Meyer-Windel (1999), which yield a S between 1.4 and 2.8 x 103 mgBB kg–1soil based on organic matter and clay contents. The sorption data and fitted isotherm are shown in Fig. 1. Most of the duplicates overlap. Since the soil was sampled from different parts of the monolith, this indicates relatively homogeneous sorption properties.
The movement of a sorbing solute is retarded relative to the mean water flow because of partitioning of the solute between the liquid and the solid phase, so that
 | [2] |
where Ct is the total concentration (mg L–1soil) and 
b the bulk density (kg L–1soil).
Sorption Kinetics
Duplicate samples of 25 g of soil and 25 mL of BB solution were used to evaluate the adsorption kinetics. Soil mixtures with two BB initial concentrations (10 and 1000 mg L–1) were shaken for 0.5, 1, 2, 4, 8, 16, and 24 h. After shaking, they were centrifuged for 15 min at 2250 g. Figure 2
illustrates that the sorbed concentration did not change with time, so adsorption was nearly instantaneous in these batch experiments (Fig. 2).
View larger version (18K):
[in this window]
[in a new window]
|
Fig. 2. Sorbed concentration dynamics for two Brilliant Blue concentrations: 10 (red squares) and 1000 mg L–1 (blue squares).
|
|
One-Dimensional Modeling
Despite its known limitations for nonlinear sorbing transport simulations (e.g., Hassanizadeh, 1996; Acharya et al., 2005), an effective one-dimensional model with constant dispersivity length was preferred to a three-dimensional model for performing the simulations, mainly because of a lack of spatial information on parameter distributions. Furthermore, when adequately parameterized, one-dimensional models have proven their utility in predicting reactive solute transport and are used for pesticide registration purposes. Main assumptions in using a one-dimensional model are that the sorption characteristics are homogeneously horizontally distributed within the monolith and that the dispersivity length adequately represents the flow field variability.
The one-dimensional CDE is used for simulating equilibrium transport in the monolith:
 | [3] |
where 
is the dispersivity length (cm), and v the average pore water velocity (cm d–1). If we assume that the sorption sites can be divided into two fractions, one fraction f with instantaneous sorption and the other one 1 – f with rate-limited sorption, then BB transport can be modeled using the following equation:
 | [4] |
where Si is the sorbed concentration on the sorption sites with instantaneous sorption and Sk on the kinetic sorption sites. Si is derived directly from the equilibrium sorption isotherm S and the fraction of instantaneous sorption sites f as
 | [5] |
and the kinetic sorption is described as a first-order rate process:
 | [6] |
where 
is the first order sorption rate coefficient (h–1).
One-Dimensional Simulations
A series of simulations of BB transport were performed with the Hydrus 1-D code (
imunek et al., 1998). The solute transport equations (Eq. [4]
–[6]) were solved using the Galerkin finite element method with a Crank–Nicholson time weighting scheme. The final solute mass balance error never exceeded 1%.
The upper boundary condition for BB transport was a 19-h-long steady-state application with constant concentration, followed by a zero-flux for 8 h. To focus on solute transport and avoid additional modeling uncertainty related to the water content and pore water velocity, the unsaturated water flow during the steady-state dye infiltration stage was not explicitly modeled but the water content profile was derived from TDR probe measurements. This was implemented in the Hydrus 1-D simulation by defining 19 soil layers with an apparent saturated hydraulic conductivity equal to the experimental flow rate Ksat = 22.7 cm d–1 and an apparent saturated water content, 
sat, that was defined based on the measured steady-state water content profile. Water contents were either interpolated between the observation depths or deduced from hydraulic considerations. At the clay layer level the water content was defined equal to the saturated water content of the sand by assuming that convergent flow induced nearly saturated conditions. Figure 3
shows the observed and "reconstructed" steady-state water content measured during chloride breakthrough at the same flow rate. Some outliers were observed in the TDR data (green open squares).
View larger version (22K):
[in this window]
[in a new window]
|
Fig. 3. Average observed water content by TDR probes under steady-state conditions (squares) and soil cores water content after drainage (crosses). Green open squares show the outliers. Continuous and dashed red lines are the simulated water content profile respectively at steady state and after drainage; continuous and dashed blue lines are the effective water content profiles at steady state and after drainage, respectively (see text for explanations).
|
|
For the 8-h drainage stage, water content and water flux time series were simulated by solving the Richards equation. Therefore, the remaining van Genuchten–Mualem hydraulic parameters (r, 
, l, n) were defined to reproduce as precisely as possible the measured water content profile after drainage, which was measured at the 104 Kopecky ring locations (crosses in Fig. 3). Unfortunately, no information was available on the dynamics of the water desaturation. This certainly adds uncertainty to the simulated water fluxes and solute redistribution during the drainage phase. However, the simulated and measured water content profiles before and at the end of the drainage phase are similar (Fig. 3), so the water redistribution, and consequently also the solute redistribution at the end of the drainage phase, can be assumed to be well represented by the simulations.
For solute transport, all the soil layers were parameterized with a dispersivity length of 3 cm (maximum observed dispersivity length in the monolith during inert BTC experiments), an average b = 1.58 g cm–3 (estimated from Kopecky rings), and the aforementioned Langmuir isotherm parameters. The largest observed dispersivity length was chosen to ensure that the total flow variability and the spreading are taken into account. Parameters and f of Eq. [5] to [7] were adapted following five scenarios summarized in Table 1. Note that apparent saturated water content profile was adapted as well for certain scenarios.
|
RESULTS AND ANALYSIS
|
|---|
BB Concentrations
The optimal model relating pixel standardized variables and BB concentrations had four variables and five parameters given with their standard in Table 2. Note that the value was found to be significant, which was not observed in other studies. On the other hand, was dismissed, and this could be due to the fact that we used artificial light to ensure the same illumination across the different depths and that the background soil color remained fairly constant with depth. The coefficient of determination of the regression model was r2 = 0.96.
Figure 4
shows the BB concentration patterns obtained from image analysis in the 12 upper monolith slices. The locations of tensiometers and TDR probes are also represented. The most remarkable feature of the BB pattern is the heterogeneity induced by a discontinuous clay layer that was observed between 0.2 and 0.28 m depth. While the concentration pattern was rather homogeneous above the clay layer, it appeared heterogeneous underneath. The patterns below 0.30 m depth are clearly an image of the shape of the discontinuous clay layer. Kung (1990) showed that coarse sands layers or densely packed fine layers could concentrate an initially unsaturated flow into irregularly spaced columns (funneling flow).
View larger version (66K):
[in this window]
[in a new window]
|
Fig. 4. Brilliant Blue lognormal concentration at 12 depths. The TDR probes are drawn in red, tensiometers in gray. Profile 1 corresponds to the bottom right-hand probes, Profile 2 to the upper probes, and Profile 3 to the left-hand bottom probes.
|
|
Since no dye tracer was observed at the outlet of the monolith, the total mass of applied BB was still in the monolith at the time of slicing and a mass balance could be performed by comparing the added tracer mass with the observed mass estimated from image analysis. The ratio between input mass and stocked dye mass was 1.07, which is good considering the numerous assumptions and simplifications made in the data processing (Forrer et al., 2000).
It is worth noting that the observed spatially heterogeneous BB pattern contrasts with the apparently local-scale homogeneous transport as assessed from inert tracer BTC experiment analysis (Javaux and Vanclooster, 2003). For the chloride BTC experiments, a longitudinal dispersivity L smaller than 1 cm was estimated at the TDR scale, while at the monolith scale L was lower or equal to 3 cm, suggesting rather homogeneous transport. Obviously, the high spatial resolution of the BB images enhances the observed transport variability. It must also be noted that the BB transport was not prone to fingering flow induced by unstable wetting front or wall flow since wetting front instability is not relevant for transport under steady-state flow conditions and no preferential flow paths were observed along the monolith walls (Fig. 4). Funneling induced by structural features appears to be the only process responsible for the flow path convergence below the discontinuous clay layer. The effect of the funneling flow was also visible in the horizontally averaged BB concentration profile which showed a pronounced tailing. In a soil with such a low dispersivity length and for a substance undergoing nonlinear sorption, which leads to a self-sharpening of the solute front, the tailing of the BB profile is a priori unexpected.
The experimental BB profile (Fig. 5
) is further characterized by two peaks at depths of 0.10 and 0.35 m located on both sides of the discontinuous clay layer and by lower concentrations near the soil surface, which may be due to the drainage after the irrigation. Kasteel et al. (2002) observed a comparable, double-peaked BB profile in their study and attributed this to flow divergence below a 10-cm-thick region where preferential flow paths were activated. Although the structural feature inducing preferential flow in our study is not a compacted plow layer, its effect on the water flow and the BB profile was similar. That is, there occurs a convergence of water into zones with a high local solute transport velocity in the layer that induces preferential flow, which accelerates the BB movement and decreases the area fraction through which transport takes place and hence decreases the averaged concentrations. Subsequently, there is a divergence of flow in the layer underneath the initiating layer, which leads to an increase in the area through which transport takes place and an increase in the horizontally averaged concentration.
View larger version (17K):
[in this window]
[in a new window]
|
Fig. 5. Observed (continuous blue line) and simulated (dashed lines) Brilliant Blue (BB) total concentration profiles considering convection–dispersion (CDE) with equilibrium (Scenario 1, in red), CDE with nonequilibrium sorption sites (Scenario 2, in green), two-site sorption (Scenario 3, in black), CDE with effective water contents and instantaneous sorption (Scenario 4, in cyan), and CDE with effective water contents and noninstantaneous, two-site sorption (Scenario 5, in magenta). The clay layer is shown in gray.
|
|
Spatial Variability of Local Concentrations
Figure 6
shows a mass fraction versus area fraction plot of the concentration distributions in the different slices. Such a plot represents the maximal fraction of dye mass that can be observed in a certain fraction of the cross-sectional area and is obtained by ranking the pixel concentrations in descending order and cumulating the concentrations and pixel areas. This is a way of characterizing the spatial solute variability in a plane with a one-dimensional graphic, analogous to the spatial solute distribution curve (SSDC) defined by De Rooij and Stagnitti (2002). A homogeneous distribution would be represented by a straight 1:1 line. It is observed that the clay layer induces important heterogeneity in the dye distribution from the 18-cm depth downward. Another interesting feature is the decrease of variability observed between 0 and 2 cm depth, which points at the effect of the upper boundary condition on the BB spatial distribution.
View larger version (30K):
[in this window]
[in a new window]
|
Fig. 6. Mass fraction of the Brillian Blue concentration distribution at different slices versus area fraction. The dashed black line represents a homogeneous distribution.
|
|
To investigate the spatial variability induced by the irrigation, we calculated the spatial variograms of the dye concentration at 0 cm depth (just below the gravel layer) for several angles relative to the grid border (Fig. 7
). Generally the plateau is reached after a 0.040-m lag distance. This is much smaller than the observation of Kasteel et al. (2005) in a silty soil field plot. The effect of the regular grid of irrigation needles on the spatial variability of local concentrations can be seen in the semivariograms. For the 90° and 0° direction the semivariance decreases a bit at the 4-, 8-, and 12-cm lag distances. In the 45° direction, this decrease can be observed at 
= 5.7 cm lag. However, the variability due to the irrigation pattern appears small compared with the total variability (
5%). Nevertheless, the redistribution occurring during the drainage phase and the soil variability may have hidden the actual variability of the solute application during leaching.
View larger version (24K):
[in this window]
[in a new window]
|
Fig. 7. Directional variograms with angles of 0° (green), 45° (red), and 90° (blue) for the upper concentration image at z = 0 cm.
|
|
One-Dimensional Modeling
Comparison between Measured and CDE-Predicted BB Profiles
The theoretical maximum total concentration estimated from the fitted sorption isotherm and the input concentration C0 = 1047 mg L–1 is, by Eq. [2], 2.11 mgBB cm–3soil. However, the maximum horizontally averaged concentration observed in the monolith is around 1 mgBB cm–3soil, just below the soil surface (Fig. 5). Moreover, as compared with a one-dimensional CDE (Eq. [3]) simulation using the batch adsorption isotherm (Eq. [1]), the shape of the observed concentration profile is far from what is expected for a nonlinear solute leaching within a uniform soil (dashed red line in Fig. 5). In addition to overestimating the maximum concentration, the simulated BB profile is delayed and has a much steeper front than observed. This difference between predicted and observed BB profiles questions the capability of models parameterized using inert tracer leaching experiments in combination with batch sorption experiments to predict the transport of a sorbing substance. It suggests that BB transport is much more influenced by the heterogeneity of the flow field than the inert tracer. It should be noted that this sequential model parameterization is typically used for practical applications.
Several processes, which were not represented in the one-dimensional CDE equation, can explain the observed BB profile: (i) flow variability induced by the macrostructure, (ii) rate-limited sorption, (iii) nonuniform distribution of the dye at the soil surface, and (iv) other microscopic processes not taken into account. In the following, these four hypotheses are reviewed in light of the available data. The aim is to investigate, by comparing model simulations with data, which of the above process hypotheses must be incorporated into the model to predict properly the BB profile in a heterogeneous subsoil. Given the need to divide the soil profile into 19 layers (with 19 parameter sets), it was impossible to perform inverse modeling, and parameters were adjusted manually.
Effect of Sorption Kinetics
As observed in Fig. 2, BB sorption in batch is almost instantaneous. However, it is known that these batch experiments may accelerate and intensify the sorption process by increasing the mixing and improving the contact between soil and solution (Schweich et al., 1983). We therefore tried to quantify the sorption rate needed to reproduce the observed BB concentration profile based on the same equilibrium adsorption isotherm.
The BB profile predicted by a CDE accounting for nonequilibrium sorption (f = 1 in Eq. [5]–[6], Scenario 2 in Table 1) with a first-order rate coefficient = 0.035 h–1 is shown by the green dashed line in Fig. 5. It should be noted that the first-order rate coefficient reflects a much slower sorption than observed in the batch experiment. Furthermore, the shape of the simulated profile does not match very well with the observed one, which probably means that adding sorption kinetics is not sufficient to explain the observed profile.
If a two-site sorption model is used (f 
1 or 0 in Eq. [5]–[6], Scenario 3 in Table 1, black dashed line in Fig. 5), the optimal parameter set is = 0.001 h–1 and f = 0.45. Since one degree of freedom is added to the one-site rate-limited sorption model, the simulation performed a bit better. The upper part of the BB profile is well reproduced, but it seems impossible to fit the lower part. Moreover, the effect of the clay layer on the simulated profile is not accounted for at all.
Inclusion of the Macrostructure in a One-Dimensional Model
As such, none of the Scenarios 1 to 3 were able to reproduce the experimental profile. This is probably due to the impact of the clay layer and the effect of the nonuniform irrigation, which may affect the flow field variability and the concentration distribution.
Nonhomogeneous irrigation actually increases the solute velocity because the water flux is focused under the irrigation needles, whereas the rest of the surface has zero flux. However, infiltration fingers were not clearly visible on the dye images, and at the 2-cm depth the effect of the irrigation variability was not visible anymore, suggesting that flow divergence quite effectively homogenized the nonuniform irrigation. The discontinuous clay layer has a similar impact on the dye distribution by inducing convergence above and divergence of the flow below the clay layer. This leads to a decrease and subsequent increase with depth of the fraction of the cross-sectional area through which transport effectively takes place.
These macroscopic features affect (i) the effective water content through which transport takes place and (ii) the effective amount of sorption sites that are accessible to the solutes. The effect of three-dimensional flow and the convergence and divergence of the flow can be represented in a one-dimensional model by adjusting the water content and the sorption isotherm with depth by a factor 
, which represents the macroscopic area through which flow effectively took place. If we consider that 50% of the surface is effectively irrigated (corresponding to an average influence diameter of 3.2 cm by needle), the effective water content represents only 50% of the measured water content. As the clay layer covers approximately 30% of the monolith cross section, also the effective water content is 30% of the total water content. To take into account the flow divergence below the soil surface, the factor was increased from 0.5 to 1 between 0 and 2 cm. The flow convergence above the clay layer is simulated by decreasing from 1 to 0.3 between 16 and 20 cm depth, and divergence by increasing from 0.3 to 1 between 20 and 46 cm depth. Therefore, for the following simulations the steady-state water content profile was changed to an effective steady state profile (z), shown in Fig. 3. As just a fraction is conducting the solute through the monolith, a fraction 1 – of the sorption sites is not accessible to the sorbing solute. Therefore, the sorption isotherm Eq. [1] must be scaled by a factor as well (this was done by adapting the parameter Cs in Eq. [1]). Transport is therefore modeled using Eq. [3] (no sorption kinetics) considering the effective water profile given in Fig. 3, and the scaled isotherm Eq. [1]. This corresponds to Scenario 4 in Table 1.
The simulated BB profile is shown in Fig. 5 (dashed cyan line). As compared with the simulation with the initial water profile, the travel depth is a bit higher, and the concentration at the soil surface lower. But, the lower concentrations in the upper 10 cm of the monolith cannot be explained by the effect of the nonideal irrigation and the macrostructure on the flow. The effect of the flow convergence and divergence due to the discontinuous clay layer at 20 cm below the monolith surface is not visible since the solute front did not reach that level. It should be noted that a simulation using an effective water content leads to an average travel time of 25 h for an inert tracer to travel to through the entire soil column. This is comparable to the observed travel time of the chloride at the same flow rate (Javaux and Vanclooster, 2003).
Combination of Macrostructure with Microscopic Effects
Considering the macrostructure and nonideal irrigation is clearly not sufficient to properly simulate the dye tracer profile. On the other hand, if only microscopic processes leading to a rate-limited sorption are considered, the solute front at the depth of the clay layer and below it cannot be described accurately, so macroscopic effects on the flow also need to be considered. Therefore, BB depth profiles were predicted by combining the macrostructure effect, modeled using an effective water content and sorption isotherm profile, with a microscopic process, that is, two-site rate-limited sorption.
Rate-limited sorption was modeled using f = 0.4 and = 0.001 h–1. The first peak and the secondary peak due to the flow divergence (macroscopic effect) are well reproduced (dashed magenta line in Fig. 5) as well as the lower than expected concentrations in the upper part of the soil monolith, which is explained by rate-limited sorption.
|
DISCUSSION
|
|---|
Including a rate-limited sorption seemed to be necessary to explain the observed BB profile, especially the lower than expected concentrations at the soil surface. However, batch sorption studies did not show evidence of a rate-limited sorption process. Other studies have already shown and explained the discrepancies existing between batch and column studies (e.g., Johnson et al., 2003), as well as the effect of the structure on nonlinear sorbing solute transport (e.g., Kasteel et al., 2002; Acharya et al., 2005). However, there are few studies on the effect of water content change on sorption rate. Maraqa et al. (1999) investigated the effect of water saturation on retardation and showed an increase of the retardation factor for lower degrees of water saturation (<70%). However, their experiment could not discriminate between lower velocity and water content diminution impacts. Gamerdinger et al. (2001) showed that rate-limited sorption was the primary mechanism at higher moisture content (>66%) in a sandy soil, and that the exchange rate was positively correlated with the flow velocity. They also showed that, at faster velocities, the sorption was lower than in batch experiments. For lower water content and velocity, the decrease of the sorption could be explained by a two-region flow, which was validated for inert tracers as well.
Estrella et al. (1993) discerned four steps in the sorption process: (i) diffusion from bulk pore water to the vicinity of mineral surface, (ii) film diffusion, (iii) intraparticle diffusion or diffusion into smaller pore sizes, and (iv) chemical sorption to the mineral surface. In our case, the chemical sorption step was quite rapid, as revealed by the batch experiment. Intraparticle diffusion at first glimpse appeared unlikely because the soil is mainly sandy and intraparticle diffusion would also have affected inert tracer transport, but this was not observed (Javaux and Vanclooster, 2003). However, it could be that BB was excluded from a smaller size pore region, even in a sandy soil, because of the large diffusion coefficient and size of BB compared with those of the inert tracer. This exclusion from smaller pores was not observed in batch experiments performed on disturbed material. The two first processes are also likely reasons for the observed rate-limited sorption in the monolith. The unsaturated conditions in the soil monolith reduce the accessibility of the sorption sites, which was obviously not the case in the batch sorption experiments. In that sense, the term chemical nonequilibrium is a bit misleading since the practical processes are "physical" (Ball et al., 1991). The aforementioned diffusive transport through fluid films from the bulk fluid toward the sorption sites on the mineral surface may be relevant for an important fraction of the less accessible sorption sites, but only represents a small fraction of the total water content, so it is not relevant for the transport of an inert tracer.
|
CONCLUSIONS
|
|---|
The BB concentration distribution obtained by image analysis after a nonsaturated leaching experiment in an undisturbed, macrostructured, subsoil monolith showed several nonideal features: (i) a deeper than expected average travel depth, (ii) an important tailing in the average BB concentration profile with a secondary peak, (iii) a lower than expected average total concentration in the upper part of the soil profile where the invading dye tracer solution was assumed to have replaced the initial solution, and (iv) a decrease of the concentration variability between the surface and the 2-cm depth. These results were in sharp contrast to the inert transport experiments that produced small dispersivity lengths and unimodal BTCs, suggesting homogeneous flow.
By comparing one-dimensional modeling results and experimental data, we investigated which processes controlled the BB transport. It was shown that:
- A two-site sorption model in combination with an effective water content and sorption isotherm that represent the effect of the three-dimensional flow induced by the macroscopic structure (i.e., a discontinuous clay layer and the regular infiltration pattern induced by the irrigation device) had to be implemented to predict the dye profile.
- The rate-limited sorption was not observed in batch sorption experiments. We postulate that this microscopic two-site sorption model was justified by the nonsaturated flow conditions, which render a fraction of the sorption sites at the grain surface less accessible for the BB solution, decreasing the sorption rate for a fraction of the sorption sites. However, more research is needed to elucidate the possible impact of water content on sorption rate in unsaturated undisturbed soil cores.
The combination of leaching experiments, statistical analyses, and one-dimensional modeling appeared to be a powerful means to point out the processes responsible for nonideal transport of a nonlinear sorbing solute in the undisturbed subsoil core. The application of the quantitative dye tracing technique was shown to give crucial information about nonideal transport processes. First, information about the structure and spatial variability of the flow process can be obtained with a high spatial resolution. Second, concentration patterns can reveal information about in situ sorption processes which cannot be observed in batch sorption experiments.
|
REFERENCES
|
|---|
- Acharya, R.C., S. Van Der Zee, and A. Leijnse. 2005. Transport modeling of nonlinearly adsorbing solutes in physically heterogeneous pore networks. Water Resour. Res. 41:W02020. doi:10.1029/2004WR003500.
- Albrecht, A., U. Schultze, P.B. Bugallo, H. Wydler, E. Frossard, and H. Fluhler. 2003. Behavior of a surface applied radionuclide and a dye tracer in structured and repacked soil monoliths. J. Environ. Radioact. 68:47–64.
- Ball, W.P., M.N. Goltz, and P.V. Roberts. 1991. Modeling the transport of solutes influenced by multiprocess nonequilibrium- comment. Water Resour. Res. 27:653–656.
- De Rooij, G.H., and F. Stagnitti. 2002. Spatial and temporal distribution of solute leaching in heterogeneous soils: Analysis and application to multisampler lysimeter data. J. Contam. Hydrol. 54:329–346.
- Estrella, M.R., M.L. Brusseau, R.S. Maier, I.L. Pepper, P.J. Wierenga, and R.M. Miller. 1993. Biodegradation, sorption, and transport of 2,4-dichlorophenoxyacetic acid in saturated and unsaturated soils. Appl. Environ. Microbiol. 59:4266–4273.
- Flury, M., and H. Flühler. 1995. Tracer characteristics of Brilliant Blue Fcf. Soil Sci. Soc. Am. J. 59:22–27.
- Flury, M., H. Flühler, W.A. Jury, and J. Leuenberger. 1994. Susceptibility of soils to preferential flow of water—A field study. Water Resour. Res. 30:1945–1954.
- Flury, M., and N.N. Wai. 2003. Dyes as tracers for vadose zone hydrology. Rev. Geophys. 41(1):1002. doi:10.1029/2001RG00010941.
- FOCUS Groundwater Scenarios Workgroup. 2000. FOCUS groundwater scenarios in the EU review of active substances. Report of the FOCUS Groundwater Scenario Workgroup. EC Document Reference SANCO/321/2000 re.2.
- Forrer, I. 1997. Solute transport in an unsaturated field soil: Vizualisation and quantification of flow patterns using image analysis. Ph.D. diss. 12476. Swiss Federal Institute of Technology, Zurich.
- 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.
- Gamerdinger, A.P., D.I. Kaplan, D.M. Wellman, and R.J. Serne. 2001. Two-region flow and rate-limited sorption of uranium (VI) during transport in an unsaturated silt loam. Water Resour. Res. 37:3147–3153.
- Hassanizadeh, S.M. 1996. On the transient non-Fickian dispersion theory. Transp. Porous Media 23:107–124.
- Haws, N.W., B.S. Das, and P.S.C. Rao. 2004. Dual-domain solute transfer and transport processes: Evaluation in batch and transport experiments. J. Contam. Hydrol. 75:257–280.
- Hinz, C. 2001. Description of sorption data with isotherm equations. Geoderma 99:225–243.
- Javaux, M., 2004. Solute transport in a heterogeneous unsaturated subsoil: Experiments and modeling. Ph.D. diss. 45. Faculté d'ingénierie biologique, agronomique et environnementale. Presses universitaires de Louvain, Louvain-la-Neuve, Belgium.
- 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.
- Javaux, M., and M. Vanclooster. 2006. Three-dimensional structure characterisation and transient flow modelling of a variably saturated heterogeneous monolith. J. Hydrol. (Amsterdam) (in press).
- Johnson, G.R., K. Gupta, D.K. Putz, Q. Hu, and M.L. Brusseau. 2003. The effect of local-scale physical heterogeneity and nonlinear, rate-limited sorption/desorption on contaminant transport in porous media. J. Contam. Hydrol. 64:35–58.
- Kasteel, R., M. Burkhardt, S. Giesa, and H. Vereecken. 2005. Characterization of field tracer transport using high-resolution images. Available at www.vadosezonejournal.org. Vadose Zone J. 4:101–111.
- Kasteel, R., H.-J. Vogel, and K. Roth. 2002. Effect of non-linear adsorption on the transport behaviour of Brilliant Blue in a field soil. Eur. J. Soil Sci. 53:231–240.
- Ketelsen, H., and S. Meyer-Windel. 1999. Adsorption of Brilliant Blue Fcf by soils. Geoderma 90:131–145.
- Kung, K.-J.S. 1990. Preferential flow in a sandy vadose zone: 1. Field observation. Geoderma 46:51–58.
- Liedl, R., and T. Ptak. 2003. Modelling of diffusion-limited retardation of contaminants in hydraulically and lithologically nonuniform media. J. Contam. Hydrol. 66:239–259.
- Mansell, R.S., S.A. Bloom, and P. Nkedikizza. 1995. Phosphorus transport in Spodosols impacted by dairy waste. Ecol. Eng. 5:281–299.
- Maraqa, M.A., R.B. Wallace, and T.C. Voice. 1999. Effect of water saturation on retardation of groundwater contaminants. J. Environ. Eng. ASCE 125:697–704.
- Olsson, J., M. Persson, J. Albergel, R. Berndtsson, P. Zante, P. Ohrstrom, and S. Nasri. 2002. Multiscaling analysis and random cascade modeling of dye infiltration. Water Resour. Res. 38(11). doi:10.1029/2001WR000880.
- Persson, M., H. Yasuda, R. Berndtsson, P. Zante, S. Nasri, and P. Ohrstrom. 2001. Modeling plot scale dye penetration by a diffusion limited aggregation (DLA) model. J. Hydrol. (Amsterdam) 250:98–105.
- Schertzer, D., and S. Lovejoy. 1987. Physical modeling and analysis of rain and clouds by anisotropic scaling multiplicative processes. J. Geophys. Res. 92:9693–9714.
- Schwartz, R.C., K.J. Mcinnes, A.S.R. Juo, and C.E. Cervantes. 1999. The vertical distribution of a dye tracer in a layered soil. Soil Sci. 164:561–573.
- Schweich, D., M. Sardin, and J.P. Gaudet. 1983. Measurement of a cation-exchange isotherm from elution curves obtained in a soil column: Preliminary results. Soil Sci. Soc. Am. J. 47:32–37.
- Seuntjens, P. 2002. Field-scale cadmium transport in a heterogeneous layered soil. Water Air Soil Pollut. 140:401–423.
- imunek, 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. IGWC-TPS70. International Ground Water Modeling Center, Golden, CO.
- Sonon, L.S., and A.P. Schwab. 2004. Transport and persistence of nitrate, atrazine, and alachlor in large intact soil columns under two levels of moisture content. Soil Sci. 169:541–553.
- Vanclooster, M., D. Mallants, J. Vanderborght, J. Diels, J. Vanorshoven, and J. Feyen. 1995. Monitoring solute transport in a multilayered sandy lysimeter using time domain reflectometry. Soil Sci. Soc. Am. J. 59:337–344.
- Vanderborght, J., P. Gahwiller, and H. Flühler. 2002. Identification of transport processes in soil cores using fluorescent tracers. Soil Sci. Soc. Am. J. 66:774–787.
- Vanderborght, J., R. Kasteel, and H. Vereecken. 2006. Upscaling of solute fluxes to the field scale using stochastic continuum transport equations. Available at www.vadosezonejournal.org. Vadose Zone J. 55 (in press).
- Young, D.F., and W.P. Ball. 1999. Two-region linear/nonlinear sorption modeling: Batch and column experiments. Environ. Toxicol. Chem. 18:1686–1693.
TOP
ABSTRACT
INTRODUCTION
