| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
a Dep. of Soil Sciences, SLU, Box 7014, 750 07 Uppsala, Sweden
b Agrosphere Institute, ICG–IV, Forschungszentrum Jülich GmbH, D–52425 Jülich, Germany
c Currently, Soil Protection Group, Insitute of Terresterial Ecology, ETH Zurich, Grabenstrasse 11a, 8952 Schlieren, Switzerland
* Corresponding author (Mats.Larsbo{at}mv.slu.se)
Received 23 September 2004.
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSREFERENCES |
|---|
Abbreviations: GLUE, generalized likelihood uncertainty estimation
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSREFERENCES |
|---|
In recent years, our knowledge of the mechanisms that generate and sustain preferential movement of water and solutes has been incorporated into several simulation models (Feyen et al., 1998; Jarvis, 1998; 
im
nek et al., 2003). Dual-permeability models divide the total soil pore space into one part (e.g., soil matrix) characterized by a large storage capacity and small flow capacity and another part (e.g., macropores) with a small storage capacity and a large flow capacity. One example of this type of model is MACRO (Jarvis, 1994), which couples classical treatments of flow and transport processes in the matrix (Richards' equation, convection–dispersion equation) to a macropore region where flow is assumed to be gravity-driven. MACRO has been widely used, both as a research tool (e.g., Larsson and Jarvis, 1999; Kätterer et al., 2001) and in management (e.g., in pesticide regulation in the EU, Forum for the Coordination of Pesticide Fate Models and Their Use, 1995), because it is physically based, numerically robust for all soil hydrological types (even for long-term simulations, i.e., decades), and is relatively parsimonious with respect to parameter requirements (imnek et al., 2003). Despite these advantages, a number of limitations of the model have also been identified since its introduction more than 10 yr ago. First, explicit numerical solutions are used for the matrix region, leading to long run times since extremely small time-steps are sometimes needed to assure stability of the numerical solutions. Second, only 22 numerical layers are allowed, restricting applications to the soil root zone, and leading to a poor spatial resolution near the soil surface. Third, the use of the Brooks–Corey retention function limits model flexibility in matching to measured hydraulic properties.
The main limitation to a much wider adoption of dual-permeability models in both research and management is the difficulty in parameterization. More parameters are required compared with classical approaches (imnek et al., 2003), some of which cannot be directly measured. The use of inverse modeling techniques may help to resolve these difficulties (Hopmans and imnek, 1999). However, the numerical limitations in MACRO noted above severely restrict the possibilities of using inverse modeling techniques to derive estimates for parameters controlling macropore flow. Roulier and Jarvis (2003a) proposed a methodology for estimating macropore flow parameters in MACRO on the basis of the global search algorithm SUFI (Abbaspour et al., 1997). However, the method is time-consuming and demands considerable computer resources, and it would not have been feasible to implement this methodology on a routine basis with the previous version of MACRO.
Here we present an improved version of the MACRO model (MACRO 5.1) (Larsbo and Jarvis, 2003) for variably saturated water flow and solute transport in structured porous media. The new features include: (i) fully implicit numerical solutions for water, heat and solute transport in the matrix, for up to 200 computational layers; (ii) use of a modified van Genuchten retention function; and (iii) inverse capabilities for parameter estimation and model calibration. First we briefly describe the model, focusing special attention on the numerical solution procedures used to couple the two flow regions. A verification of the numerical solutions is then presented, in which simulation results are compared with analytical solutions for water flow in macropores (Germann, 1985). Our objectives were (i) to test the ability of the model to simulate high time-resolution measurements of water flow and nonreactive (Cl–) solute transport (including effluent and resident concentrations) in transient microlysimeter experiments and (ii) to investigate the identifiability of four critical model parameters regulating macropore flow using the GLUE procedure (Beven and Binley, 1992).
|
MATERIALS AND METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSREFERENCES |
|---|
Soil Water Flow
Richards' equation is used to calculate vertical water fluxes in the matrix:
 | [1] |
/
(m–1) is the differential water capacity, (m3 m–3) is the volumetric water content, (m) is the soil water pressure head, t (s) is time, z (m) is depth, K (m s–1) is the unsaturated hydraulic conductivity, and Si (s–1) are source–sink terms accounting for water exchange with macropores, drainage, and root water uptake.
The use of Eq. [1] to calculate water flows in the macropore domain is problematic; a major reason for this is the lack of information concerning () close to saturation. Therefore capillarity is assumed to be negligible in the macropores, so water flow is driven by gravity only (i.e., /z = 0). The governing equation for water flow in macropores is the kinematic wave equation (Germann, 1985):
 | [2] |
ma (m3 m–3) and Kma (m s–1) are the macropore water content and hydraulic conductivity, respectively.
Hydraulic Properties
In macroporous soils, hydraulic conductivity increases very rapidly across a small pressure head range as saturation is approached (Clothier and Smettem, 1990; Jarvis and Messing, 1995). In MACRO, a "cut and join" approach is used to define the matrix–macropore hydraulic functions (Jarvis et al., 1991). Thus, a user-defined pressure head, b (m), partitions the total porosity into matrix pores and macropores, while a corresponding water content, b (m3 m–3), and hydraulic conductivity, Kb (m s–1), represent the saturated state of the soil matrix.
Soil water retention in the matrix is calculated using a modified form of van Genuchten's (1980) equation (Vogel et al., 2001; Fig. 1)
:
 | [3] |
vg (m–1) are shape parameters (where mvg is set equal to 1 – 1/nvg); mi (m3 m–3) is the micropore water content; r (m3 m–3) is the residual water content; and s* 
is a "fictitious" saturated water content, obtained by extrapolating the fitted water retention function to zero pressure head. It should be noted that s* does not represent the actual saturated water content in the model, which is separately defined by the user to reflect macroporosity.
|
 | [4] |
b) is the effective water content at micropore saturation given by replacing with b in Eq. [3].
The hydraulic conductivity function in the macropores is given as a simple power law of the macropore degree of saturation, Sma:
 | [5] |
Water Uptake by the Matrix
In the absence of gravity, Richards' equation can be recast as a diffusion equation, with gradients in water content as the driving force. In MACRO, imbibition of water from macropores into an unsaturated matrix is treated as a first-order approximation to this water diffusion equation, assuming a rectangular-slab geometry for the aggregates (van Genuchten and Dalton, 1986; Booltink et al., 1993):
 | [6] |
w is a scaling factor introduced to match the approximate and exact solutions to the diffusion problem (Gerke and van Genuchten, 1993). The scaling factor w varies with the initial water content and hydraulic properties, but not strongly (Gerke and van Genuchten, 1993; Jarvis, 1994), so for simplicity w is set within the program to an average value (Jarvis, 1994). The effective water diffusivity is given by
 | [7] |
b (m2 s–1) and Dmi (m2 s–1) are the water diffusivities at the saturated matrix water content and the current matrix water content, respectively, and where Sma is introduced to account for incomplete wetted contact area between the two pore domains. In the Mualem–van Genuchten model, Dmi is given by (van Genuchten, 1980)
 | [8] |
is a "fictitious" saturated hydraulic conductivity obtained by extrapolating the matrix conductivity function to zero pressure head. Db is given by setting S in Eq. [8] to Smi(b). It can be noted that the true saturated hydraulic conductivity of the soil (see Eq. [5]) will usually be larger than Ks* as a result of the effects of macropores.
Equation [6] is treated as a source term to Eq. [1] and a sink term to Eq. [2]. Water transfer in the reverse direction (from the matrix to macropores) occurs instantaneously when the pressure head in the matrix exceeds b, following the fundamental physical principle governing the filling of pores when the water-entry pressure is exceeded. The way this is resolved is further described below in the section Generation of Macropore Flow.
Solute Transport
Solute transport is calculated using the convection–dispersion equation with source–sink terms, Ui (kg m–3 s–1), to represent a wide range of processes, including mass exchange between flow domains, kinetic sorption, solute uptake by plants, biodegradation, and lateral leaching losses to drains and/or regional groundwater:
 | [9] |
mi(m) (m3 m–3) is the accessible water content accounting for anion exclusion, q (m s–1) is the water flow rate, and D (m2 s–1) is the dispersion coefficient calculated as the sum of an effective diffusion coefficient and a dispersion term. In the macropores, an equivalent approach is used to calculate transport, except that dispersion is assumed to be zero and only equilibrium sorption is considered. Equilibrium sorption partitioning is calculated using the Freundlich isotherm.
Solute Uptake by the Matrix
The source–sink term for mass transfer of solute between matrix and macropores, Ue (kg m–3 s–1), is given by a combination of a diffusion component and a mass flow component:
 | [10] |
Numerical Solutions
At each time step, vertical water and solute fluxes are first calculated in the matrix. Updated values of water storages are used to determine the excess amount of water routed to the macropores. Water fluxes in the macropores are then calculated, and the solute concentrations are derived by an implicit elimination technique.
Equation [1] is solved using the iterative procedure proposed by Celia et al. (1990). Details of its implementation into MACRO are given in Larsbo and Jarvis (2003). Here, we focus on aspects of the numerical solution that are particularly relevant to dual-permeability models, especially the surface boundary condition and the coupling to the macropore region. We also briefly describe the numerical method used to solve Eq. [2] describing flow in macropores.
Boundary Conditions
Water will flow into surface-vented macropores if the rainfall intensity exceeds the infiltration capacity of the matrix. The net rainfall, R (m s–1), at the soil surface is partitioned into soil evaporation, Es (m s–1), an amount infiltrating into the matrix, Imi (m s–1), and an excess amount flowing into the macropores, Ima (m s–1), depending on the infiltration capacity of the matrix, Imax (m s–1), calculated from Darcy's Law. Thus, if R – Es < Imax, then the top boundary condition for the matrix is defined by a known flux (Neumann condition) equal to the net rainfall minus soil evaporation. Otherwise, the soil surface boundary condition is given by a known tension (Dirichlet condition).
The actual soil evaporation rate, Es, is limited either by the potential evaporation rate or the maximum possible rate of supply of water to the soil surface (Feddes et al., 1974). The latter is calculated explicitly by iteratively substituting Mualem–van Genuchten's model into Darcy's Law and evaluating the flux as a function of the unknown pressure head at the surface and the known pressure head at the first node.
Generation of Macropore Flow
As shown in Fig. 1, a pressure head, b, and the corresponding water content, b define the division between flow domains. If the water content in the matrix exceeds b, the excess water will drain into the macropores. However, the macropores only have a finite storage capacity for water, Cma (m3 m–3), equal to
 | [11] |
max [m3 m–3], max [m], see Fig. 1) defines the maximum amount of water allowed in the matrix at each timestep (max = b + Cma). When solving Richards' equation, the pressure heads are allowed to increase above max only if the macropores are also saturated. In this situation, the differential water capacity in the matrix (Eq. [1]) is set to zero for pressures above max. For pressure heads above b, the hydraulic conductivity is set to Kb. Once Richards' equation has been solved, all water exceeding b is instantly removed from the matrix and added to the macropores.
Water Flow in Macropores
The local water balance in the macropores (Eq. [2]) is solved implicitly for Sma using a bisection method ("interval halving") for each layer, in turn, from the soil surface downward. If the flow capacity of both matrix and macropores is exceeded in any layer, then the local water balance cannot be solved (oversaturation develops in the macropores). In this case, the excess water is added to the macropore storage in the layer(s) above, and the water fluxes between layers are corrected accordingly.
Convection–Dispersion Equation
The convection–dispersion equation is solved using a Crank–Nicholson finite difference scheme. An iterative, fully upstream weighted procedure is used which minimizes overshoot problems and oscillations that may sometimes be encountered in the presence of steep concentration gradients. However, upstream weighting introduces considerable numerical dispersion, which is minimized by an empirical correction:
 | [12] |
The surface boundary condition for the matrix is a solute flux, Jmi (kg m–2 s–1), which can be either positive or negative depending on the mass balance:
 | [13] |
is the solute concentration in water routed into the macropores, calculated by assuming instantaneous local equilibrium and complete mixing of incoming water and solute with the water and solute stored in a shallow surface soil layer or "mixing depth" (Steenhuis and Walter, 1980; Jarvis, 1994). The surface flux into the macropores is simply given by the second term on the right-hand side of Eq. [13].
Model Verification
The model has been compared with analytical solutions of one-dimensional water flow and solute transport problems to verify the accuracy of the numerical solutions (Vanderborght et al., 2005). The model fit to the analytical solutions for a range of flow and transport problems was generally very good. Here we present a comparison between simulations and analytical solutions to the kinematic wave equation (Eq. [2] and [5]) for macropore flow (Germann, 1985). A 2-m-deep soil profile divided into 200 numerical layers was simulated. Figure 2
shows that the numerical solution is in good agreement with the analytical solution, even though some numerical dispersion of the wetting and draining fronts is apparent. It can be noted that such an accurate numerical solution to the macropore flow problem could, paradoxically, result in a poorer match to observations, in the case of the wetting front advance, and the initial stage of drainage. This is because, in reality, flow takes place in a range of macropores of different sizes and continuity, so that some dispersion of the wetting and draining fronts is usually found (Germann, 1985).
|
s* = 0.502 m3 m–3, nvg = 1.13, and vg = 0.63 m–1. Setting b to –10 cm gave b = 0.500 m3 m–3. The microlysimeters were set up in the laboratory allowing free drainage at the base through two layers of wire mesh (hole size = 1 mm2) and one layer of cloth to prevent particle migration (hole size 0.4 mm2). Irrigation was supplied through pump-driven atomizer spray devices located 1.2 m above the center of each microlysimeter. The atomizers deliver a fine spray mist of low kinetic energy, so the original structure of the soil surface is retained throughout the experiment.
Four irrigations of 17.1, 11.8, 18.6, and 18.8 mm of filtered natural rainwater were supplied to the microlysimeters, the first lasting 3 h, and the following three each lasting 2.5 h. The first two irrigations were performed with a 4-d interval between them. The third irrigation commenced approximately 7 d after the second, and the final irrigation started 36 h later. Before the third and fourth irrigations, a KCl solution was applied to the soil surface using a small hand-held sprayer, which supplied doses of 74.8 and 81.8 g Cl m–2 in the two irrigations, respectively, in 1.42 mm of water. Samples of water outflow were collected at time resolutions varying between 3 and 20 min depending on the water outflow rates. Thirty-six hours after the fourth irrigation, the soil was excavated from the microlysimeters in five layers, each 3 cm thick. To extract Cl–, a known quantity of distilled water was added to each sample, followed by end-over-end shaking for 24 h, and filtration to remove any particulate matter. Chloride concentrations in soil extracts and water samples were measured by flow injection analysis. Only one of the four microlysimeters was chosen for the modeling exercise since all four showed similar patterns of percolation and Cl– transport.
The GLUE procedure (Beven and Binley, 1992) was applied to the microlysimeter experiment to evaluate the possibilities to identify four key model parameters determining the strength of preferential flow: the saturated matrix hydraulic conductivity (Kb), the macroporosity (ma,s [m3 m–3]) calculated as the difference between s and b, the kinematic exponent (n*), and the effective diffusion pathlength (d) (Table 1). These parameters were chosen either because they are impossible to measure directly or because they were thought to be among the most sensitive based on prior experience with the model (Dubus and Brown, 2002). The GLUE methodology explicitly recognizes the underlying limitations of environmental models by accepting that all models and measurements are to some extent wrong. Consequently, we cannot expect to find one unique optimal parameter set for a specific model application. Many parameter sets may equally well represent the observations according to some measure of goodness-of-fit. The GLUE procedure starts off with predefined parameter distributions and uncertainty limits. The parameter uncertainty limits (Table 1) were based on a preliminary analysis performed on much wider uncertainty intervals (Kb, ma,s, and d) and prior experience with the model (n*). Sampling should be made from parameter distributions that are consistent with expected model responses to changes in parameter values (Beven and Binley, 1992). We transformed d to d–2 since mass exchange between the flow domains is inversely proportional to the square of the diffusion pathlength (Eq. [6] and [10]). Kb was sampled from log-transformed values based on a preliminary analysis of model response. Uniform distributions were used for the remaining parameters. The remaining parameters in the model were either set to known (i.e., diffusion coefficient for Cl–, 1.9 x 10–9 m2 s–1) or assumed values (e.g., matrix dispersivity 5 cm, mixing depth 1 mm), or estimated from pedotransfer functions included in MACRO 5.1 (saturated hydraulic conductivity = 31.9 mm h–1). The initial water content and a constant potential evaporation rate were determined through preliminary calibrations. A "lysimeter" bottom boundary condition (zero potential for saturated conditions, no inflow for unsaturated conditions) was used in accordance with the experimental setup.
|
We ran 20000 simulations generated by Latin hypercube sampling from within the parameter uncertainty limits and calculated EF values for all simulations. Simulations with overall model efficiencies larger than a predefined threshold EFtot value were considered "acceptable." Even though GLUE was designed to identify acceptable parameter sets, information on individual parameters can be obtained from cumulative distributions of model efficiencies for acceptable simulations (Beven and Freer, 2001). These distributions give information on the degree of parameter conditioning. Conditioned parameters, as opposed to fitted parameters, result in model outputs that "honor" all or parts of the measured data (Abbaspour et al., 1999). Those parameters with cumulative distributions that differ the most from the initial uniform distributions have been most conditioned by the process. A steep "S"-shaped cumulative distribution means that most acceptable simulations have parameter values in a limited part of the initial uncertainty interval. Hence, the cumulative distributions give information on the possibilities to reduce the uncertainty in a parameter by some calibration method. The choice of threshold value may strongly influence the cumulative distributions. We calculated two sets of cumulative distributions for threshold values of EFtot > 0.5 and EFtot > 0.3. The threshold value was subtracted from each EFtot value before the calculation of cumulative distributions. Uncertainty ranges defined by the minimum and maximum values of the simulated output variables from simulations with EFtot larger than the threshold were compared with measured data. Clearly, these uncertainty ranges depend on the choice of goal function and on the threshold value.
|
RESULTS AND DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSREFERENCES |
|---|
|
|
|
The EFtot values from the GLUE simulations are presented in Fig. 6
. Parameter sets with EFtot values >0.5 are concentrated in a limited part of the parameter space for d–2, whereas, for log(Kb) and ma,s they are scattered over larger parts of the parameter space and for n* they are scattered over the whole parameter space. It can be noted that values of d–2 larger than approximately 0.05 mm–2 (d < 5 mm) resulted in small EFtot values compared with the best GLUE simulation. This means that a certain degree of preferential flow is needed to simulate the experiments well.
|
ma,s, and d–2 or to reduce their uncertainty intervals by inverse techniques, while identification of n* may not be possible. It is not surprising that we encountered problems in simultaneously identifying both ma,s and n* since the water flow through the macropores is dependent on both these parameters (Eq. [5]). A negative correlation between ma,s and n* was also found for the simulations with EFtot larger than 0.5 (Table 2). Using response surface analysis of d and n*, Roulier and Jarvis (2003b) showed that n* could not be identified in their similar microlysimeter experiments because of lack of sensitivity. Table 2 also shows a negative correlation between log(Kb) and d–2. This can be explained by the fact that the effect of high macropore infiltration [small log(Kb) values] on Cl– leaching, can be compensated by a fast exchange between the pore domains (large d–2 values).
|
|
ma,s, and n*. However, the advantages of microlysimeters (e.g., relating parameter values to horizon properties) would then be lost. Alternatively, ma,s could be estimated directly from specific yield measurements to allow better identification of n*, while Kb could be estimated through tension infiltrometer measurements in the field to allow better identification of d. These parameters could also be retained in the calibration procedure, with uncertainty domains and distributions constrained by the measurements. Flow interruption experiments have been used successfully for determination of the degree of preferential flow in mobile–immobile models (Reedy et al., 1996). This type of experiment allows complete equilibration between the flow domains and may further improve the conditioning of d. |
CONCLUSIONS |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSREFERENCES |
|---|
Data from a transient microlysimeter tracer breakthrough experiment (high resolution measurements of percolation and effluent concentrations and a final resident concentration profile) indicating strong preferential flow were accurately reproduced by model simulations (EFtot = 0.62). A GLUE analysis for four parameters regulating the degree of preferential flow showed that this type of experiment contains sufficient information for conditioning of the matrix hydraulic conductivity, the macroporosity, and the diffusion pathlength, indicating that these parameters would be identifiable in inverse modeling approaches based on microlysimeter experiments. However, the kinematic exponent was not satisfactorily conditioned by the data, which was attributed primarily to correlation with the macroporosity. The computationally intensive GLUE analysis was feasible with the new version of the model because of significantly reduced run times compared with previous versions.
|
ACKNOWLEDGMENTS |
|---|
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTS AND DISCUSSIONCONCLUSIONSREFERENCES |
|---|
imnek. 1999. Review of inverse estimation of soil hydraulic properties. p. 643–659. In M.Th. van Genuchten et al (ed.) Characterization and measurement of the hydraulic properties of unsaturated porous media. Proc. Int. Workshop. Univ. of California, Riverside.imnek, J., N.J. Jarvis, M.Th. van Genuchten, and A. Gärdenäs. 2003. Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone. J. Hydrol. (Amsterdam) 272:14–35.This article has been cited by other articles:
|
|
 |
|
  N. Jarvis Near-Saturated Hydraulic Properties of Macroporous Soils Vadose Zone J., November 26, 2008; 7(4): 1256 - 1264. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
S. Panday and P. S. Huyakorn MODFLOW SURFACT: A State-of-the-Art Use of Vadose Zone Flow and Transport Equations and Numerical Techniques for Environmental Evaluations Vadose Zone J., May 27, 2008; 7(2): 610 - 631. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 M. Larsbo, K. Fenner, K. Stoob, M. Burkhardt, K. Abbaspour, and C. Stamm Simulating Sulfadimidine Transport in Surface Runoff and Soil at the Microplot and Field Scale J. Environ. Qual., May 1, 2008; 37(3): 788 - 797. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
T. H. Skaggs, N. J. Jarvis, E. M. Pontedeiro, M. Th. van Genuchten, and R. M. Cotta Analytical Advection Dispersion Model for Transport and Plant Uptake of Contaminants in the Root Zone Vadose Zone J., November 20, 2007; 6(4): 890 - 898. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
H. H. Gerke, J. Dusek, T. Vogel, and J. M. Kohne Two-Dimensional Dual-Permeability Analyses of a Bromide Tracer Experiment on a Tile-Drained Field Vadose Zone J., August 23, 2007; 6(3): 651 - 667. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
A. Boivin, J. Simunek, M. Schiavon, and M. Th. van Genuchten Comparison of Pesticide Transport Processes in Three Tile-Drained Field Soils Using HYDRUS-2D Vadose Zone J., June 21, 2006; 5(3): 838 - 849. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||
