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


     

Published online 27 April 2006
Published in Vadose Zone J 5:515-528 (2006)
DOI: 10.2136/vzj2005.0056
© 2006 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 (3)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Javaux, M.
Right arrow Articles by Vanclooster, M.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Javaux, M.
Right arrow Articles by Vanclooster, M.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Javaux, M.
Right arrow Articles by Vanclooster, M.
Related Collections
Right arrow Solute Transport Models
Right arrow Dispersion
Right arrow Vadose Zone Processes and Chemical Transport

ORIGINAL RESEARCH

Three-Dimensional Modeling of the Scale- and Flow Rate-Dependency of Dispersion in a Heterogeneous Unsaturated Sandy Monolith

M. Javauxa,*, J. Vanderborghta, R. Kasteela and M. Vancloosterb

a Agrosphere, ICG-IV, Forschungszentrum Jülich, D-52425 Jülich, Germany
b Dep. of Environmental sciences and Land Use Planning, Université Catholique de Louvain, Croix du Sud, 2 Bte 2, B-1348 Louvain-la-Neuve, Belgium
* Corresponding author (m.javaux{at}fz-juelich.de)

Received 18 April 2005.



   ABSTRACT
TOP
 EXECUTIVE SUMMARY
 ABSTRACT
INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Spatial heterogeneity in soil hydraulic properties strongly affects solute transport. However, the level of precision needed in the description of spatial variability to properly reproduce the solute mixing regime is still an open question. We investigated this problem by analyzing observed inert solute transport using three-dimensional simulations with different levels of complexity in the material description. The scale- and flow-dependency of the dispersivity was first characterized from a series of leaching experiments during unsaturated steady-state flow in a heterogeneous sandy monolith. The structure and hydraulic property variability within the monolith was investigated as well by means of an exhaustive survey of the monolith, and by intensive soil core sampling allowing for hydraulic characterization. In this study, three three-dimensional models are constructed, involving several levels of complexity. In Case I, only the macrostructure variability is represented. In Case II, scaling factors encoding the spatial variability in the hydraulic properties of the sandy matrix are implemented. In Case III, an anisotropy factor for hydraulic conductivity is added to the macrostructure and the microheterogeneity of the sand matrix. Results show that microheterogeneity is needed to reproduce qualitatively the scale- and flow rate-dependency of the transport parameters. Despite the elaborate effort devoted to the structure characterization, no model was fully capable of reproducing observed solute transport in the monolith and at the outlet.

Abbreviations: BTC, breakthrough curve • EC, electrical conductivity • TDR, time domain reflectometry


   INTRODUCTION
TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
AQUANTITATIVE UNDERSTANDING of the impact of the spatial heterogeneity of a porous medium on solute transport is key to developing models applicable to large-scale soil and groundwater pollution problems. When a discrete hierarchical description of heterogeneity (Cushman et al., 2002) is used, the relevant description of the heterogeneity is a function of the scale of observation and of modeling. In the following we define microheterogeneity as the heterogeneity that cannot be described in a deterministic way because the scale over which the parameters vary is too small compared with the large-scale modeling domain. Since the distribution of parameters at each location cannot be defined in a deterministic sense, a stochastic representation of the heterogeneity is necessary. The spatial structure of this type of heterogeneity can be represented, for example, by a spatial covariance function. By comparison, macroheterogeneity refers to structures of a larger spatial extent, which can be parameterized in a deterministic way within the modeling domain.

The impact of spatial heterogeneity on macrodispersivity previously has been studied numerically for media with different levels of complexity. First, some authors investigated the influence of microheterogeneity on flow and transport in partially saturated isotropic Miller-similar, two-dimensional, macroscopically homogeneous media (Birkholzer and Tsang, 1997; Roth, 1995; Roth and Hammel, 1996). The main conclusions from these studies are that (i) microheterogeneity promotes transport through a complex pattern of channels, (ii) the hydraulic structure is flow dependent, and (iii) a critical flow state exists where flow heterogeneity (and longitudinal dispersivity) is minimal. However, approximate analytical solutions of the stochastic flow and transport equations indicate that the existence of this critical state depends on the correlation between the different hydraulic parameters (Russo, 1998). The critical flow rate exists in Miller-similar media or in media where scaling factors for hydraulic conductivity, fK, and pressure head, fh, are negatively correlated. For noncorrelated or positively correlated fK and fh, flow becomes increasingly more heterogeneous when the saturation degree decreases.

Cirpka and Kitanidis (2002) increased the medium complexity by accounting for cross-correlation between the heterogeneous fields of unsaturated parameters in a three-dimensional loamy soil. In their case, the Ksat and {alpha} parameters of the Mualem–van Genuchten functions were positively correlated (corresponding to a negative correlation between fK and fh). They showed that for small flow rates, local-scale dispersion and local-scale mixing become relatively more important, so that the solutes are effectively mixed between different regions. Given that the minimal variability in the flow field is at a lower infiltration rate for loamy soil than for a sand, transverse diffusion has more time to smooth out concentration differences, which causes the macroscopic dispersivity to increase with flow rate and depth.

Several other authors (Abdou and Flury, 2004; Khaleel et al., 2002; Tseng and Jury, 1994; Ursino et al., 2000) analyzed the effect of medium anisotropy on transport. In a deterministic framework, anisotropy refers to a given hydraulic parameter being different in different directions due to pore-scale structures. In a stochastic framework, anisotropy is related to the spatial structure of an hydraulic parameter, which shows a larger correlation range in one direction than in the other.

Ursino et al. (2000) extended the concept of Miller-similar media to deterministically anisotropic microheterogeneous structures and included a change in the anisotropy factor with saturation. They characterized a switching point at which anisotropic media behave like an isotropic column. The effect of the stochastic anisotropy was investigated by Khaleel et al. (2002), who generated a two-dimensional medium with independent hydraulic parameters characterized by correlation lengths identical for all parameters, but different in the horizontal and vertical directions. They found an increase in the longitudinal macrodispersivity when the soil becomes drier. All these studies compared flow and transport parallel and perpendicular to bedding, in which case the soil tended to behave as a parallel column system or a layered system, respectively (Ursino et al., 2000). Khaleel et al. (2002) also reported a larger dispersivity and mixing length for flow parallel to bedding.

Working on a two-dimensional artificial anisotropic heterogeneous sand tank consisting of three soil textures arranged in small quadratic layers, Ursino et al. (2001a, 2001b) introduced macroheterogeneity and showed that low saturation levels produced considerable heterogeneity in the flow pattern and preferential flow. They stressed the fact that saturation-dependent macroscopic anisotropy is an essential element for modeling transport in unsaturated media. Subsequently, Ursino and Gimmi (2004) tried to model their observations after having characterized the structure pattern by image analysis and knowing the characteristic hydraulic functions of each soil type. Unfortunately, they were not able to simulate appropriately the experimental results. This was attributed to small-scale microscopic transport processes, which are difficult to consider explicitly in the frame of a macroscopic continuum approach.

Hammel et al. (1999) used a hierarchical description of heterogeneity, with a coarse component defined by the deterministic pattern of soil horizons, and a local fine-scale component described with a correlated random field of scaling factors conditioned on measured hydraulic functions within the horizons. They compared simulated and observed bromide distributions in a heterogeneous field soil and showed that only simulations with stochastic local and deterministic large-scale hydraulic variations gave appropriate results. Deurer et al. (2001) used a comparable modeling concept for a layered sandy soil. They defined four reference parameter sets corresponding to the four horizons and used Vogel scaling factors for characterizing the microheterogeneity of each layer. Although their model reproduced a similar channel-like solute distribution, quantitative agreement with the data was poor; in particular, they could not replicate the preferential flow of solutes in the deeper horizons.

Hence, when a hierarchical description of heterogeneity is used, macrodispersion will depend on the way in which micro- and macroheterogeneity interact and are described within the transport model. With respect to microheterogeneity, previous studies showed that macrodispersion will depend on (i) the local-scale spatial distribution of flow and transport properties, (ii) the correlation between hydraulic parameters or scaling factors, and (iii) the degree and type of anisotropy incorporated in the model. Following Vogel and Roth (2003), for inert solute transport, this microheterogeneity can be dealt with by means of lumped effective parameters within a macroheterogeneous medium. For saturated media, several studies have proposed a rigorous framework for upscaling heterogeneities to the larger scale (Attinger et al., 2004; Rubin et al., 1999). Yet, it is presently still not clear how macro- and microheterogeneity interact in unsaturated media and how such media should be described in efforts to model macroscopic inert solute transport. In addition, few experimental studies have been presented for testing the performance of the different methods for implementing heterogeneity in a hierarchical framework as a mean to model macroscopic inert solute transport in real undisturbed soils.

We attempt to evaluate different approaches for implementing heterogeneity in a hierarchical heterogeneous medium by comparing inert solute transport processes assessed experimentally in an unsaturated heterogeneous monolith using three-dimensional simulations with different levels of complexity in the material description. In a recent study, Javaux and Vanclooster (2003) performed a comprehensive study to characterize the scale- and flow-dependent solute transport processes in a 0.5-m3 unsaturated heterogeneous monolith. In addition, an exhaustive survey of the structure of the medium allowed them to assess the hydraulic property variability (Javaux and Vanclooster, 2006). The objective of our study was to test the possibility of predicting effective solute transport and its scale- and flow rate-dependency based on the observed (micro- and macro-) structure and hydraulic properties of the heterogeneous medium.


   MATERIALS AND METHODS
TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Physical Model
The experimental study was conducted using a cylindrical monolith of 0.5 m3 (1-m height, 0.77-m i.d.) extracted from a sand quarry, 10 m below the normal ground surface. The studied subsoil is a Tertiary marine deposit, so-called Bruxellian sand, covering the central part of Belgium.

In the laboratory, the monolith was equipped with 12 time domain reflectometry probes (TDR) horizontally inserted along three vertical transects, spaced 120° apart, with two of the transects having three probes (located at the 0.15-, 0.45-, and 0.75-m depths) and the third transect having six probes (located at the 0.15-, 0.30-, 0.45-, 0.60-, 0.75-, and 0.90-m depths). Four tensiometers were located at depths of 0.1, 0.3, 0.6, and 0.95 m. Flow and electrical conductivity (EC) meters allowed continuous measurement of outflow and the EC at the monolith outlet. All sensors except the tensiometers recorded data automatically.

Solute Transport Experiments
Eleven steady-state chloride breakthrough curve (BTC) experiments involving Dirac pulses were performed with the monolith at different flow rates, ranging from 2 to 100% of the saturated conductivity: eight under unsaturated conditions and three under saturated conditions, with a zero-tension bottom boundary condition. Various features of the experiments are given in Javaux and Vanclooster (2003). The main findings were that the apparent dispersivity and the mixing regime was dependent mostly on the scale of observation and far less on the flow rate. Subsequently to the chloride tests, a dye tracer experiment was performed with Brilliant Blue, after which the monolith was sliced to observe the macroscopic spatial heterogeneity (Javaux, 2004).

Assessment of the Structure
The heterogeneity of the subsoil core was experimentally characterized at two spatial levels. At the macroscopic scale, we delimited three texturally contrasted bodies in the monolith—a discontinuous clayey layer and a stony layer embedded in a sandy matrix. The bodies were identified by means of cross-section image analysis and by observation during the monolith destruction.

The microstructure of the sandy matrix was characterized by an exhaustive sampling of the monolith. We randomly extracted 104 undisturbed soil cores of 100 cm3 at seven depths within the monolith in the sand body. Based on multistep outflow experiments and suction–water content pairing, the hydraulic properties of the samples at these 104 locations were obtained by inverse modeling of the one-dimensional Richards equation. The resulting data permitted an assessment of the stochastic textural variability of the sand. Details about the experiments and results are given in Javaux and Vanclooster (2006).

Soil Macrostructure
Figure 1 gives the location and extent of the 9-cm-thick clay layer, located at the 20-cm depth, and the discontinuous stony layer embedded at the 73-cm depth (2–3 cm thick). The three soil types were considered to be adequately characterized by the van Genuchten–Mualem hydraulic functions. The parameter sets of the clay and stone layers were estimated from the limited data we had since no direct information was available about their hydraulic behavior. For example, the saturated water content and hydraulic conductivity of the clay layer were estimated with the Hypres pedotransfer function (Wösten et al., 1998), based on textural analysis of the clay layer (60% clay, 40% silt, and 0.6% organic matter). Parameters of the stone layer were defined to keep the stone layer saturated and with a relatively low hydraulic conductivity (much lower than the sand matrix) because of cementation. The parameters of the sand material corresponded to the reference parameters obtained by the scaling procedure. The reference hydraulic parameters for the different soil types are reported in Table 1.


Figure 1
View larger version (40K):
[in this window]
[in a new window]
  		Fig. 1. Macrostructure (red: stones, yellow: clay) in the monolith and position of the TDR probes. TDR of the first transect (cyan), second transect (green), and third transect (blue).

 

View this table:
[in this window]
[in a new window]
  		Table 1. Soil hydraulic parameters.

 
Sand Local Variability Parameterization
Javaux and Vanclooster (2006) assessed the sand microvariability using the Vogel-similarity concept (Vogel et al., 1991). The variability was encoded in the following four scaling factors which relate the parameters of van Genuchten–Mualem reference hydraulic curves (noted hereafter with *) to the characteristic parameters at other locations x:

Formula 1[1a]

Formula 2[1b]

Formula 3[1c]

Formula 4[1d]
where f{theta}r, fK, fh, and fSe are the scaling parameters for, respectively, the residual water content ({theta}r), the saturated conductivity Ksat, the suction and the degree of saturation and where {theta}sat is the saturated water content, and {alpha} is a parameter of the van Genuchten hydraulic function inversely related to the inflection point of the water retention curve.

The four scaling factors described were estimated at 100 locations within the sand body in the monolith (four outliers were detected and dismissed for further analysis), for which experimental cross- and autovariograms were obtained after subtraction of a depth trend. Subsequently, isotropic cross- and autovariograms (Fig. 2 ) were parameterized using a nested Gaussian-exponential model defined as

Formula 5[2]
where aij1 and aij2 are the optimized sills for models g1 and g2, respectively. The exponential model g1 is defined for short range (b1 = 0.25 m) as

Formula 6[3]
and where g2 is the gaussian model for longer range (b2 = 0.9 m):

Formula 7[4]
We used a linear model of coregionalization to parameterize the ten auto- and cross-semivariograms. More details are given in Javaux and Vanclooster (2006).


Figure 2
View larger version (22K):
[in this window]
[in a new window]
  		Fig. 2. Experimental and modeled auto- and cross-variograms of the four scaling factors.

 
Three conditional fields of parameters were generated using the sequential Gaussian simulation algorithm to simulate the heterogeneous spatial distribution of the four scaling factors within the monolith. However, only one simulation was used for simulating the nine steady-state leaching experiments. The two other simulations were used to assess uncertainty related to simulations at one given flow rate. Conditional simulations incorporated the known scaling factors at 100 locations to simulate possible outcomes at unsampled locations. To avoid excessive computing time, the simulations were performed on a grid with cells of 0.03-m height and 0.017-m width (smaller than 5 times the apparent correlation length). Finally, these simulations were linearly interpolated to match the finer three-dimensional model grid (see next section). All processing was done in MATLAB 6.0 (Mathworks Inc., Natick, MA) using the BMElib toolbox (Christakos et al., 2002) for the statistical part. Since textural variability was only assessed for the sandy matrix, scaling factors in the stony and clay layers were set equal to 1 (no variability).

Figure 3 shows a series of cross sections of the simulated scaling factors within the monolith. Note that the macrostructure is not presented in the graphic for the sake of clarity. Given the large correlation range of most factors, hot spots can be observed; these hot spots are related to local observed hydraulic parameters. Note also that inverse correlation between fK and fh, or between fSe and f{theta}r for instance, are readily apparent.


Figure 3
View larger version (62K):
[in this window]
[in a new window]
  		Fig. 3. Horizontal cross sections in the monolith showing distribution of the four scaling factors.

 
Anisotropy
Because of the sedimentation process, the subsoil showed some slight layering. The fact that no anisotropy appeared in the directional semivariograms of the scaling factors (not shown here) is due to the sampling height of the soil cores (5 cm), which contained more than one layer. However, some anisotropy in the local conductivity tensor was expected. To test this hypothesis, a factor of 10 between the horizontal and vertical conductivity will be assumed in several scenarios (to be discussed more below).

We decided not to incorporate a moisture-dependent anisotropy factor mainly for two practical reasons: (i) no information was available on this dependency and (ii) it would lead to much larger computation times. However, this type of anisotropy could be easily implemented in further simulations (e.g., Raats et al., 2004).

Three-Dimensional Model
We assumed that, at the continuum scale, water flow in soil is adequately described by the Richards equation:

Formula 8[5]
with h the pressure head (L), t the time (T), {theta} the water content (L3 L–3), K the conductivity tensor (L T–1), and x1 the vertical coordinate. We similarly assume that, at the same scale, the convection dispersion equation is also valid

Formula 9[6]
with C1r the total resident concentration (M L–3), v (L T–1) and D (L2 T–1) the velocity and dispersion tensors. Components of the dispersion tensor are defined as

Formula 10[7]
where {delta}ij is the Kronecker delta function, {tau} is a tortuosity factor, Dd is the molecular diffusion coefficient (L2 T–1), and {lambda}l and {lambda}t are the local longitudinal and transverse dispersivities (L), respectively. In our study, we neglected molecular diffusion (Dd = 0), while the local transport parameters are considered constant for the whole porous medium: {lambda}l = 5 x 10–3 m and {lambda}t = 5 x 10–4 m. The value for the longitudinal dispersivity was chosen in accordance with the lowest dispersivity observed at the TDR scale at the upper level of the monolith (Javaux and Vanclooster, 2003). We also assumed that the effective transport volume for solute corresponded with the total local water content.

The numerical model SWMS_3D (Simunek et al., 1995) was used to simulate three-dimensional flow and transport in a rigid isotropic porous medium using variable flux boundary conditions at the top boundary. This model solves the Richards equation for saturated–unsaturated water flow using the finite element method with linear elements. To eliminate undesired oscillation in the numerical solution due to convection-dominated transport, we used upstream weighting in the solution of the transport equation and required that the product of the Courant by Peclet numbers be ≤2.

A MATLAB program was written to automatically generate three-dimensional mesh grids for a cylindrical flow domain compatible with the SWMS_3D code. The soil was discretized into 471 765 hexahedral elements (automatically subdivided in tetrahedral by the program) of 0.009975-m width and 0.01-m height (i.e., 105 slices of 4493 elements). This fine discretization was chosen to avoid numerical dispersion. A preliminary test was performed by simulating several leaching experiments for a homogeneous soil with the same discretization and a constant dispersivity of 5 x 10–3 m. The dispersivity obtained by fitting the outlet BTC always corresponded to that same 5 x 10–3 m.

We note that original code considered three scaling factors, namely fK, fSe, and fh, to represent the variability of the hydraulic functions. Minor modifications were made to implement f{theta}r as the fourth scaling parameter.

Flow and Transport Boundary Conditions
For the transport simulation we assigned an initial concentration of zero throughout the flow domain. The bottom boundary conditions were zero-tension for flow and a zero- diffusive flux for transport. Table 2 summarizes the upper flow boundary conditions for the nine simulations. Nine different steady-state flow rates, jw0, including one saturated flux, were considered. These rates corresponded to the 11 (eight under nonsaturated conditions and three replicates for saturated conditions) applied flow rates for the BTC experiments (Javaux and Vanclooster, 2003). The leaching experiment consisted of homogeneously releasing the solute during a short time period at the top face of the column during steady-state flow. The duration of the pulse was different for each flux to match the various experimental conditions (Table 2).


View this table:
[in this window]
[in a new window]
  		Table 2. Boundary conditions for the numerical steady-state solute leaching experiments. The three last columns give the calculated averaged water content at steady-state for the three simulated cases.

 
Simulation Scenarios
To compare the effects of microvariability and macrostructure on solute transport, three cases were evaluated. In Case I, only the macrostructure was modeled, without any microvariability. In Case II, the macrostructure and microstructure (textural variability of the sand) were taken into account. Case III groups both micro- and macrostructure effects with a certain degree of anisotropy for the hydraulic conductivity (factor 10 between the horizontal and vertical hydraulic conductivity), to take into account the effect of the slight layering resulting from the sedimentation process.

Effective Models and Virtual TDR Probes
To estimate an apparent dispersivity for the entire monolith for each flow rate, we fitted the simulated outlet BTCs with a one-dimensional analytical solution of the convection–dispersion equation. This equation was used to give apparent parameters that we could use for comparison with the experimental results by Javaux and Vanclooster (2003) using the same analytical solutions. The dispersivity obtained from the outlet BTC will be referred to further as the apparent macrodispersivity {lambda}La.

Similarly to the methodology used by Javaux and Vanclooster (2003), we also followed the concentrations at nodes corresponding to the location of TDR probes in the real monolith. Between 200 and 250 nodes covering the TDR sampling window were averaged for each time step to simulate an apparent BTC at the position of the TDR probes (Fig. 1). This was done for 12 locations described above in the Physical Model section. Note that the probes themselves and their influence on flow lines were not simulated in the three-dimensional model. To analyze the simulated TDR BTCs, we transformed the concentrations into time-normalized values Crt*(z,x,t) (Vanderborght et al., 1996).

Effective dispersivity and velocity estimates at each TDR location were obtained by fitting the one-dimensional convection–dispersion equation to the local data. The apparent dispersivity estimated at the scale of TDR probes will be referred as the TDR-scale longitudinal dispersivity {lambda}LTDR. At the three depths where three TDR probes were inserted together, we also calculated BTCs by averaging the pixels corresponding to the three TDR probes, which gave a better approximation of the macrodispersivity at these depths. It should be noted that this fitting procedure with one-dimensional CDE solutions of the was used to obtain apparent parameters characterizing the shape of the observed BTCs, but that we did not assume that the actual effective transport process in the monolith was mostly convective–dispersive as such.


   RESULTS AND DISCUSSION
TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
CONCLUSIONS
 REFERENCES
 
Steady-State Flow Simulation
Impact of Heterogeneity on Steady-State Flow
In Fig. 4 , we compare monolith cross sections of the water content ({theta}), pressure head (h), and relative magnitude of the flow vector [jw(x)/jw0) at steady state with jw0/Ksat = 0.036 (Exp. 2) for the three cases. Streamlines were also computed from the three-dimensional velocity vector. The impact of the macrostructure on these variables is clear, irrespective of the level of complexity that was implemented. Convergent–divergent flow paths are found in the vicinity of the discontinuous clay layer. Using the scaling factors (Cases II and III) made the water content, the pressure head, and the flow fields more heterogeneous and the streamlines more tortuous. However, although preferred flow pathways can be identified, one may not be able to define clear channel structures, probably because of the large correlation length of fK (Fig. 2) relative to the monolith diameter. When local anisotropy is added, streamlines become even more tortuous due to the higher lateral conductivity. This also tends to homogenize the pressure heads at the same depth. The effect on the spatial distribution of the water content was relatively small.


Figure 4
View larger version (56K):
[in this window]
[in a new window]
  		Fig. 4. Comparison of the pressure head (first column), water content (second column), logarithm of the local flow velocity and streamlines (third column) under steady state with jw0 = 2.19 x 10–7 (Exp. 2) for the three cases: Case I, only macrostructure (first row); Case II, macrostructure plus sand hydraulic parameter variability (second column); and Case III, macrostructure plus sand hydraulic variability and anisotropy for hydraulic conductivity (third row). Streamlines (in blue) were computed from linear interpolations of the velocity field.

 
Comparison with Observations
Javaux and Vanclooster (2006) investigated the capabilities of the three-dimensional models to reproduce transient flow behavior of the monolith. They showed that the outflow and pressure head time series were simulated reasonably well with and without considering microheterogeneity. To compare the three scenarios with our data, we used steady-state pressure head profiles for the nine flow rates. Figure 5 shows the average pressure head profiles and the range of pressure heads at a given depth. The pressure head variability was less when an anisotropic local hydraulic conductivity was assumed than for Case II. Notice also that the discrepancy between the models with (Cases II and III) and without (Case I) microheterogeneity tended to increase with increasing flow rate.


Figure 5
View larger version (42K):
[in this window]
[in a new window]
  		Fig. 5. Pressure head profiles at steady state for the nine flow rates (Table 2): observed data (open circles) and average pressure head profiles for Cases I (red line), II (green line), and III (blue line). Color-similar filled areas give the extent of the pressure head variability.

 
Comparison between data and simulations revealed some inconsistencies for the lowest and highest flow rates in particular. These inconsistencies can be attributed to problems with the parameterization of the hydraulic reference curve parameters {alpha} and n of the sand and of the conductivity curve of the clay layer (Javaux and Vanclooster, 2006). Ranking of the performance of the three models is further complicated by the fact that they performed relatively well for certain experiments but failed for others. Overall, however, Cases II and III appeared to be slightly better.

Overview of Solute Transport in Heterogeneous Media
Cross sections of the normalized resident concentrations of Exp. 5 (Table 2) are shown in Fig. 6 . The characteristic features described hereafter for this experiment were very much the same for all nine experiments. The effect of the three types of heterogeneity implemented in the monolith on solute plume distortion is clearly visible in the figures. For Case I, the plume shape is affected by macrostructure only, while no smaller-scale channeling is visible. The channeling effect in the Vogel-similar medium (Case II) is apparent for all flow rates on both the transversal and longitudinal cross sections. Adding anisotropy for the conductivity (Case III) leads to more redistribution horizontally such that the concentrations in the "bypassed" regions are larger and the concentration contrast in a horizontal cross section is smaller for the anisotropic case than for the isotropic case.


Figure 6
View larger version (47K):
[in this window]
[in a new window]
  		Fig. 6. Vertical and horizontal cross sections of relative solute concentration distributions in the monolith for a flow rate jw0 of 1.64 x 10–6 m s–1 (Exp. 5) after an irrigation dose of 13.16 cm. Upper right: depth-averaged relative concentration for Case I (red line), Case II (green line), and Case III (blue line). The dotted white lines show the cutting planes.

 
The normalized resident concentration distributions, calculated by summing the relative concentrations at each depth, are given in the same figure. Dispersion generated without microheterogeneity (Case I) is obviously less, although Cases II and III appear very similar, so that the effect of anisotropy on the one-dimensional averaged resident concentration profile is less obvious.

Note also that a secondary peak appears above the clay layer, as was also observed in the Brilliant Blue profile obtained by image analysis (Javaux, 2004). This secondary peaks are due to the small velocities and the divergence of flow above the clay layer. The resulting delays in solute transport induce tailing in the BTCs. On the other hand, a decrease in concentration can be observed at the depth of the stone layer. This reflects the acceleration in the solute particles produced by flow convergence in the sand between the stones.

To visually check the uncertainty related to the stochastic simulations of the van Genuchten parameters of the sand body, we conducted two other simulations and simulated flow and transport for the flow rate corresponding to Exp. 7 (Table 2). Figure 7 shows the simulations with and without anisotropy for the three simulated correlated random fields. Although three simulations are not enough to estimate any uncertainty distribution, the results show that the conditional simulations led to very similar transport process, especially for the isotropic medium. When anisotropy was added, the discrepancies were larger. This is due to the large horizontal conductivity, which allows the solute particles to bypass the less conductive soil parts more easily.


Figure 7
View larger version (50K):
[in this window]
[in a new window]
  		Fig. 7. Vertical cross sections of relative solute concentration distribution in the monolith for a flow rate jw0 of 2.63 x 10–6 m s–1 (Exp. 7) after an irrigation dose of 10.51 cm. Each column represents one different stochastic simulation. First line is for the isotropic media (Case II) and the second line for the anisotropic medium (Case III).

 
Scale- and Rate-Dependency of Simulated and Observed TDR-Scale Parameters
Apparent Velocity Profiles
Figure 8 shows the simulated and observed apparent velocity profiles obtained from the TDR probes located on Transect 1 (Fig. 1) for several flow rates. The measured apparent velocity profiles are somewhat curved by higher values below the clay layer and a decrease toward the inlet and outlet of the soil column. This curvature in the distributions was well reproduced by the three models. Cases II and III are again quite similar, while deviating from Case I. Since the microstructure contains much more uncertainty (e.g., local hydraulic conductivities are only known at the local sampling points and are simulated between sampling points using conditional Gaussian simulation and linear interpolation), deviations between the simulated and measured velocity profiles are mostly a result of this uncertainty.


Figure 8
View larger version (27K):
[in this window]
[in a new window]
  		Fig. 8. Apparent velocity obtained from TDR-probe BTCs of the first profile versus depth for Case I (red lines), Case II (green lines), Case III (blue lines), and observations (black lines) for five flow rates.

 
Apparent Dispersivity Profiles
Figure 9 compares the depth-dependency of the TDR-scale {lambda}LTDR obtained from simulations with the experimentally derived {lambda}LTDR. The error bars indicate the average and standard deviation of {lambda}LTDR at a given depth for different flow rates.


Figure 9
View larger version (23K):
[in this window]
[in a new window]
  		Fig. 9. The TDR-scale {lambda}L–depth relationships for Case I (top, in red), Case II (middle, in green), and Case III (bottom, in blue). Observed {lambda}L distributions are given in black. Probes from Profile 1 are represented by open squares, Profile 2 with open circles, and Profile 3 with open triangles (see Fig. 1 for a description of probe locations). Error bars show the standard deviation due to flow rate.

 
The observed {lambda}LTDR obtained from the three TDR probes located at the 0.15-m depth were almost independent of the flow rate, varying between 0.6 to 1 cm. This is well reproduced by the three models. The values should be compared with the constant local dispersivity {lambda}l = 0.5 cm, which was implemented in the model. Below the clay layer at around the 0.2-m depth, the observed {lambda}LTDR increased up to 1 cm. This is again well simulated with the three models, even without accounting for microheterogeneity. The variability and higher flow dependency generated by the convergence–divergence flow paths around the clay layer is better observed at the 0.45-m depth, where three TDR probes sampled the monolith cross section. Notice from Fig. 9 that Case II produced larger differences in {lambda}LTDR between TDR locations. When local anisotropy was added for the conductivity (Case III), this variability decreased (i.e., {lambda}LTDR was quite similar among the three probes at the 45-cm depth).

The three models of heterogeneity represented the observed trends very well: an increase in {lambda}LTDR between 0 and 0.30 m depth, stabilization around {lambda}LTDR = 1 cm between 0.30 and 0.90 m, and a small decrease at the 0.90-m depth. However, Case III performed better in that some overlap occurred between the observed and predicted {lambda}L distributions for all TDR probe locations. Generally, averaged {lambda}LTDR was overestimated by Cases II and III and underestimated by Case I. On the other hand, the rate-dependency of {lambda}LTDR (expressed by the extent of the error bar) was much lower for Case I, whereas Cases II and III were more sensitive to the flow rate, with this sensitivity being very similar to the observed one.

Mean Behavior on the Monolith Scale: Scale- and Rate-Dependency
Apparent Macrodispersivity at the Outlet: Impact of the Flow Rate
Figure 10 shows the apparent macrodispersivity values {lambda}La estimated from the simulated outlet BTC for the nine different flow rates, and for the three models of heterogeneity. The experimental monolith-scale dispersivity was overestimated with Cases II and III and underestimated with Case I heterogeneity. The general shape for Cases II and III, however, matched the experimental results much better.


Figure 10
View larger version (21K):
[in this window]
[in a new window]
  		Fig. 10. Effective longitudinal dispersivity estimated from simulated outlet breakthrough curves for Case I (red), Case II (green), Case III (blue), and experiments (dots).

 
The overestimation of the effective dispersivity at the outlet contrasts with the relatively good approximation of the local dispersivity within the monolith (see Fig. 9 and also the next section). It is possible that the correlation length of the conductivity in the vertical direction was overestimated, especially when assuming an isotropic microstructure. If the real structure was anisotropic, with a larger correlation range in the horizontal direction, then this would lead to an overestimation of the correlation scale in the vertical direction, and therefore, of the longitudinal dispersivity. Another explanation is that {lambda}l may be overestimated and {lambda}t underestimated. These local parameters were considered constant in our study, although some studies reported a dependency on the saturation degree (Haga et al., 1999; Padilla et al., 1999; Sato et al., 2003). On the other hand, the experimental bottom boundary conditions of the monolith may have had some effect on the measured BTC by faciltating additional mixing. Moreover, numerical dispersion could be another reason for {lambda}La overestimation.

We found almost no impact of the flow rate on the simulated {lambda}La for Case I, except a slight decrease at the larger flow rates. Cases II and III produced comparable shapes, with minimum {lambda}La values at a critical flow rate corresponding to Exp. 5 (jw0/Ksat sand = 0.27) for Cases II and III. These results involving a minimum dispersivity at a critical flow rate or suction head were expected for microheterogeneous media characterized by a negative correlation between fh and fK (Roth, 1995), as in our case (Fig. 2). However, the bottom flow boundary conditions (seepage face) and the macrostructure implemented in the models did not allow the water content to change much with jw0. It is relatively surprising that {lambda}La was larger for Case III than for Case II. Indeed, Ursino et al. (2000) investigated the effect of anisotropy on transport in a Miller-similar medium and noticed a slightly smaller longitudinal dispersivity when anisotropy was added. However, the difference between {lambda}La for Cases II and III is relatively negligible compared with the uncertainty associated with the fitting procedure. Furthermore, given the length of the monolith and the extent of the macrostructure and of correlation lengths of the microstructures, the asymptotic steady state for the transport was clearly not reached in our study.

Figure 11 shows the relative velocity cross sections for several flow rates. Due to the bottom boundary condition, the lower parts of the sections did not vary much with flow rate. At the clay layer, the higher velocity region resulting from convergence of streamlines remains for all flow rates. The changes in the relative velocity patterns are not really evident, except in the upper part of the monolith, with a switch between more and less conductive regions beyond the critical velocity as expected for Miller-similar media (Roth, 1995).


Figure 11
View larger version (81K):
[in this window]
[in a new window]
  		Fig. 11. Steady-state log-relative velocity cross sections for five different flow rates.

 
Effective Dispersivity at the Monolith Scale: Impact of the Travel Depth
To evaluate monolith-scale solute transport, we averaged the three BTCs obtained with the TDR probes at three depths (0.15, 0.45, and 0.75 m; see Fig. 1) to obtain "monolith-scale" BTCs, from which apparent macrodispersivities {lambda}La were derived. This was done for both the simulations and the experimental results. Figure 12 compares, for three flow rates, the depth-dependency of the monolith-scale macrodispersivity {lambda}La. Notice that the dispersivity obtained from the outlet BTC at the 1.05-m depth was also added in the figure.


Figure 12
View larger version (15K):
[in this window]
[in a new window]
  		Fig. 12. Monolith-scale dispersivity profile. Crosses or circles represent the average of the dispersivity obtained at a given depth for all flow rates; error bars give the standard deviation. Data, Case I, Case II, and Case III are in gray, red, green and blue, respectively.

 
Javaux and Vanclooster (2003) previously found a scale effect between the TDR-scale and monolith-scale dispersivities and a change in the apparent mixing regime between scales. Cases II and III reproduced this scale effect as well, with larger {lambda}La than {lambda}LTDR values. The depth-dependency of {lambda}La is also observed here for Cases II and III, but not for Case I. Indeed, since no correlated random field of conductivity existed for Case I, variations in {lambda}La with depth are directly related to the macrostructure only, and no scale effect is expected. The discrepancy between the experimental and Case II {lambda}La values increased with depth. Adding anisotropy (Case III) tended to decrease {lambda}La within the monolith, which produced a better match with the data. However, the outflow {lambda}La was again overestimated for Cases II and III, as explained above.

Retardation at the Outlet
Figure 13 shows the mechanical retardation factor, Rma, defined as the ratio between the velocity of water (vw = jw0/<{theta}>) and the apparent velocity of solute obtained by BTC fitting. The experimental data generally indicate a slight acceleration in the solute velocity (Rma < 1). This is also predicted with all models. As noted, Case I predicted a much faster leaching at lower flow rate, while at higher flow rates, Case II produced a higher leaching velocity. At the lower flow rate, it was the macrostructure that accelerates transport by restricting the effective transport volume. At the larger flow rates, the high-velocity regions tended to transport relatively large solute mass and thus accelerated the solute transport. Notice, however, that the Rma for all models tended to 1 with increasing flow rates.


Figure 13
View larger version (22K):
[in this window]
[in a new window]
  		Fig. 13. Mechanical retardation factor at the outlet for Case I (red circles), Case II (green circles), Case III (blue circles), and experiments (+).

 

   CONCLUSIONS
TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
REFERENCES
 
In this study we tried to reproduce BTCs that were measured during leaching experiments in a sandy monolith, on the basis of an exhaustive description of the spatial distribution of hydraulic parameters. The spatial distribution was represented using a deterministic macrostructure, a stochastic microstructure conditioned on local measurements, and by considering local anisotropy of the hydraulic conductivity. The objective was to properly predict the scale- and rate-dependency of the dispersivity and the velocity. Furthermore, we tried to elucidate the relative contribution of microstructure, macrostructure, and anisotropy on the variability of the transport of a chloride plume in the sandy monolith at several flow rates.

For this case study, the macrostructure alone could not explain the scale- and flow rate-dependency of the transport process. Inclusion of the stochastic microstructure allowed us to simulate the transport process in a qualitative sense. The agreement could be further improved by including local anisotropy in the hydraulic conductivity (Fig. 9), which results from microstructures in the medium (i.e., thin layers associated with sedimentation). Indeed, given that the layering was smaller than the scale of the soil cores used for characterization of the hydraulic functions, anisotropy could not have been observed by geostatistical analysis. This is the main reason we implemented only microscopic anisotropy in the conductivity tensor of the three-dimensional model rather than using anisotropic correlation lengths for generating the scaling factor spatial field.

Scenarios with microstructure revealed the existence of a critical flow rate where the dispersivity was minimum. This is in line with numerical simulations for media with negative correlation between fK and fh (Birkholzer and Tsang, 1997; Roth, 1995; Roth and Hammel, 1996). However, both numerical and experimental results showed that the seepage face bottom boundary conditions combined with the macrostructure tended to decrease the effect of the flow rate on the transport in a medium with a microstructure. This probably explains why the critical flow rate was closer to Ksat than in previous studies.

Despite the enormous efforts devoted to the characterization of macro- and microheterogeneity, and the controlled boundary conditions, we were not able to construct a three-dimensional model that fully reproduced the monolith hydraulic behavior. We found that substructure characterization is needed to predict the scale- and flow-dependency of the dispersivity. Still, despite the good performance of the Case III model in predicting the scale- and rate-dependency of the dispersivity inside the monolith, a discrepancy remained in the predictions at the outlet. The accurate determination of the reference hydraulic functions, the cross correlation and statistical parameters that characterize the microstructure heterogeneity, and the way in which anisotropy is implemented in the model still remain important sources of uncertainty and could explain deviations between the simulated and measured transport parameters.

As compared with larger scale transport processes, the effect of the macrostructure in the monolith experiments was likely reduced by the imposed particular lateral, surface, and bottom boundaries. These boundary conditions hinder lateral redistribution of flow that would have occurred at larger scale by these structures. Analysis of observed concentrations of the large-scale experiment suggests that these structures at the larger scale may indeed play an important role (Javaux and Vanclooster, 2004a, 2004b).

More fundamentally, this study raises the question about which level of prediction can be reached with numerical models, and whether or not real-world complexity can ever be captured in enough detail. More research is needed on the comparison between models and experiments for porous media with known structure and well-characterized hydraulic behavior. Furthermore, it would probably be worthwhile to define up to which level a model should adhere to reality, depending on the scale of the simulation, the imposed boundary conditions, and medium complexity. Finally, research of this type can give us better insights into the optimal balance between more detailed description of the porous media structures involved and model performance.


   REFERENCES
TOP
 EXECUTIVE SUMMARY
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 




This article has been cited by other articles:


Home page
 Vadose Zone JHome page
H. Hardelauf, M. Javaux, M. Herbst, S. Gottschalk, R. Kasteel, J. Vanderborght, and H. Vereecken
PARSWMS: A Parallelized Model for Simulating Three-Dimensional Water Flow and Solute Transport in Variably Saturated Soils
Vadose Zone J., April 9, 2007; 6(2): 255 - 259.
[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 (3)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Javaux, M.
Right arrow Articles by Vanclooster, M.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Javaux, M.
Right arrow Articles by Vanclooster, M.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Javaux, M.
Right arrow Articles by Vanclooster, M.
Related Collections
Right arrow Solute Transport Models
Right arrow Dispersion
Right arrow Vadose Zone Processes and Chemical Transport

HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
The SCI Journals Agronomy Journal Crop Science
Journal of Natural Resources
and Life Sciences Education		 Soil Science Society of America Journal
Journal of Plant Registrations Journal of
Environmental Quality		 The Plant Genome