Published online 16 December 2005
Published in Vadose Zone J 5:59-76 (2005)
DOI: 10.2136/vzj2005.0008
© 2005 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
ORIGINAL RESEARCH
Inverse Dual-Permeability Modeling of Preferential Water Flow in a Soil Column and Implications for Field-Scale Solute Transport
J. Maximilian Köhnea,c,*,
Binayak P. Mohantya and
Jirka 
im
nekb
a Dep. of Biological & Agricultural Engineering, Texas A&M Univ., Scoates Hall, College Station, TX 77843-2117
b Univ. of California, Dep. of Environmental Sciences, 900 University Av., A135 Bourns Hall, Riverside, CA 92521
c currently, Institute for Land Use, Univ. Rostock, Justus-von-Liebig Weg 6, 18059 Rostock, Germany
* Corresponding author (max.koehne{at}uni-rostock.de)
Received 14 January 2005.
 |
ABSTRACT
|
|---|
The question of whether or not soil hydraulic parameters of dual-permeability models (DPM) can be properly identified by inverse analysis of preferential water flow data has not been resolved to date. We applied a DPM based on two coupled Richards' equations to compare the performance of inverse and forward simulations of laboratory preferential flow data. Infiltration and drainage experiments were conducted using a repacked loam soil column (80 cm long, 24-cm diameter) containing a cylindrical sand region (2.4-cm diameter) as the preferential flow path (PFP) along its central axis. The forward DPM water flow simulations relied on hydraulic parameters for the matrix and the PFP as determined by means of separate infiltration and drainage experiments on loam and sand columns, respectively. One inverse DPM approach relied on standard (lumped) observations of infiltration and outflow, while the other included outflow through the matrix and the PFP. Both inverse approaches provided accurate matches of bulk infiltration and outflow, but the outflow out of the matrix and the PFP could only be described when fitting the DPM to region-specific outflow data. The practical implication of this finding for predicting solute transport in natural soils was evaluated. An observed tile-drainage hydrograph was used for inverse hydraulic DPM parameter estimation, followed by fitting the solute transfer coefficient and the dispersivity for simulating Br– tracer concentrations. This sequential fitting procedure was successful for hydrograph simulation but unsuccessful for Br– breakthrough simulation. Simultaneous hydraulic and transport parameter estimation considerably improved the approximation of Br– concentrations. This study shows that a hydrograph alone is not sufficient for inverse identification of soil hydraulic DPM parameters. Simultaneously employing hydrograph and solute breakthrough data may facilitate identification of hydraulic and transport DPM parameters to characterize preferential flow and solute transport.
Abbreviations: DPM, dual-permeability model • IC, initial conditions • LBC, lower boundary conditions • PFP, preferential flow path • SWCC, soil water characteristic curve • TDR, time domain reflectometry
|
INTRODUCTION
|
|---|
DUAL-PERMEABILITY MODELS are increasingly used for analysis of preferential solute transport in structured soils that contain PFPs, such as fractures or macropores, both on the laboratory column scale (e.g., Gwo et al., 1995, 1996; Allaire et al., 2002a, 2002b; Castiglione et al., 2003; Greco, 2002; Kätterer et al., 2001) and on the plot or field scale (e.g., Jarvis et al., 1991; Andreu et al., 1994; Larsson and Jarvis, 1999a, 1999b; Kohler et al., 2001). A problem receiving considerable attention recently is the determination of the many DPM parameters required for simulating preferential water flow and solute transport. While the soil hydraulic functions for DPM can be independently estimated from water retention and hydraulic conductivity measurements with double-hump curve shapes (e.g., Mohanty et al., 1997; Köhne et al., 2002), this approach cannot provide the parameters of water transfer between matrix and macropore domains that are also required. Inverse simulation is an alternative way to estimate DPM parameters that are otherwise difficult to obtain. However, questions about the uniqueness of such optimized parameters still exist (imnek et al., 2001). Several investigations on preferential water flow mechanisms in soils with PFPs were reported (e.g., Booltink and Bouma, 1991; Tofteng et al., 2002; Logsdon, 1995), but only few preferential water flow experiments have been analyzed in detail for evaluation of the region-specific hydraulic parameters related to the matrix or the PFPs (e.g., imnek et al., 2001; Castiglione et al., 2003; Köhne and Mohanty, 2005). While most studies employed a mixture of measurements, estimates, and manual calibration of DPM parameters to simulate preferential solute transport (e.g., Jarvis et al., 1991; Jabro et al., 1994; Andreu et al., 1994; Ahuja et al., 1995; Larsson and Jarvis, 1999a, 1999b; Malone et al., 2001; Gerke and Köhne, 2004), inverse parameter optimization schemes have been employed as well (e.g., Schwartz et al., 2000; Abbaspour et al., 2001; Kätterer et al., 2001).
Traditionally, the parameters of (nonpreferential) flow and transport models are estimated sequentially. First the soil hydraulic parameters are estimated using hydraulic observations, followed by the estimation of transport parameters using solute concentration data (Jacques et al., 2002; Abbasi et al., 2003). However, comparisons between this sequential fitting procedure and simultaneous estimation of hydraulic and transport parameters from combined data of water flow and solute concentrations rarely showed the sequential scheme to be superior (e.g., Abbasi et al., 2003), while other studies found the simultaneous scheme to be more efficient (Medina and Carrera, 1996; Mishra and Parker, 1989; Inoue et al., 2000). For example, imnek et al. (2002) used experimental data of Inoue et al. (2000) and showed that although the final parameter estimates for both sequential and simultaneous optimizations were very similar, the confidence intervals for the simultaneous optimization of soil hydraulic and solute transport parameters, as well as the final value of the objective function, were smaller. It remains unclear whether the sequential parameter estimation scheme can be employed for the inverse simulation of preferential water flow and solute transport in the vadose zone, since water flow observations may not always contain sufficient information for proper inverse identification of hydraulic DPM parameters (imnek et al., 2001). A recent study showed that despite an excellent fit of the tile-drainage hydrograph, a two-dimensional dual-porosity (mobile–immobile) model may fail to describe observed tracer breakthrough (Haws et al., 2005).
These findings from previous studies suggest that (i) hydrographs may not carry enough information to distinguish between local water flow paths with distinctly different water flow velocities and (ii) while domain-specific flow velocities may not be critical to describe water flow, they may become significant for simulating preferential solute transport. Although reasonable, the above conclusions are only partly supported by experimental data, since these data were obtained rather indirectly at the field scale, where many additional confounding unknowns are present, such as soil heterogeneity, off-site contributions, catchment boundaries, and initial and boundary conditions. Hence, a closer examination of the underlying causes of the above findings seems warranted.
The primary objective of this study was to investigate if DPM hydraulic parameters for matrix and PFP regions can be identified based on standard soil hydraulic observations, such as infiltration, outflow, pressure heads, and water contents of a well-controlled column experiment, or if additional domain-specific outflow through matrix and PFP must be included. The secondary objective was to evaluate the practical implications of the DPM hydraulic parameter identification for field-scale preferential solute transport prediction. Observations of drainage flux and preferential Br– tracer breakthrough in a tile-drained field between 1991 and 1994 (Lennartz et al., 1999) were simulated, and we compared the performances of sequential and simultaneous inverse estimation of soil hydraulic and solute transport DPM parameters.
|
MATERIALS AND METHODS
|
|---|
Physical Soil Properties
The particle size distributions of the sandy loam soil (henceforth referred to as loam) and coarse and medium sand are listed in Table 1. The loam, sampled from arable land near College Station, TX, was air-dried and ground. Additionally, 30 smaller soil columns (7.5 cm in diameter and height) were packed with loam and the medium and coarse sand (10 columns per material) to measure the soil water characteristic curves (SWCC) between saturation and –500 cm pressure head according to Klute (1986). Based on the 10 sets of data for each material, the average water retention functions were determined by fitting soil hydraulic functions (van Genuchten, 1980) using the RETC model (van Genuchten et al., 1991). Furthermore, the saturated hydraulic conductivity, Ks, was measured on the same soil columns (Klute and Dirksen, 1986). The results for Ks and the SWCC are provided in Table 2.
View this table:
[in this window]
[in a new window]
|
Table 1. Particle size distribution, organic carbon content (OC), and bulk density ( ) of the soil materials used in the homogeneous soil columns (i.e., loam, medium-sand, and coarse-sand columns) and in the columns with preferential flow paths (PFP).
|
|
View this table:
[in this window]
[in a new window]
|
Table 2. Soil hydraulic van Genuchten's (1980) parameters and the range of 95% confidence limits obtained by inverse solution of the Richards' equation (imnek et al., 1998) based on experimental infiltration and drainage data for homogeneous loam and medium-sand and coarse-sand columns. The parameters were used in the forward dual-permeability simulations. Corresponding parameters are also given for the soil water characteristic curve (SWCC) measured at static equilibrium and for Ks at steady-state flow (SSF).
|
|
Flow Experiments with the Homogeneous Soil Columns
Laboratory soil columns 16 cm in diameter and 40 cm high were packed with loam, coarse sand, and medium sand, respectively. These air-dry soil materials were packed using a piston compactor at 3-cm increments to obtain relatively homogeneous bulk densities between 1.6 and 1.7 g cm–3. The bulk densities theoretically correspond to porosities ranging between 0.35 and 0.4, if a particle density of 2.65 g cm–3 is assumed. After slowly saturating the columns with deionized water from the bottom to the top to avoid air inclusions, outflow and infiltration experiments were conducted to determine the soil hydraulic parameters. A tension infiltrometer with a 16-cm-diameter infiltration disc (i.e., matching the soil column) was used for infiltration at zero pressure head. At the column bottom the soil was in contact with a nylon membrane supported by a densely perforated PVC plate. The lower (seepage) boundary was subject to adjustable suctions controlled by a vacuum system with an automated outflow monitoring system. Pressure heads were recorded across the soil column with tensiometers at 10-, 20-, and 30-cm depths, while water contents at the 10- and 30-cm depths were measured using time domain reflectometry (TDR) probes of 8-cm length (the same probes were used in the preferential flow experiments described below).
Experiments were conducted for the following initial and lower boundary conditions (IC and LBC): (i) infiltration into the loam column (IC: –140 cm, LBC: 0 cm), (ii) three-step outflow out of the initially saturated loam column (LBC: 0, –60, and –120 cm), (iii) infiltration into the medium-sand column (IC: –40 cm (top) to –20 cm (bottom), LBC: –15 and 0 cm), (iv) drainage from the initially saturated coarse-sand column (LBC: 0, –30, and –60 cm). Infiltration into coarse sand was measured as well; however, unstable (fingering) flow was observed, and the results were excluded from further analysis.
Preferential Flow Experiments
Our soil column setup for studying preferential flow was described in detail in a recent study (Köhne and Mohanty, 2005). Only those setup features are mentioned here that are relevant to the present study, or that were modified. The setup consisted of a large loam soil column (24.4-cm diameter, 80-cm length, henceforth referred to as matrix) with a cylindrical sand-filled PFP (2.4-cm diameter) at the center. A 5-cm-high cylindrical flow separator at the bottom was used to prevent water exchange between the matrix and the PFP due to boundary effects. The PFP region was implemented by vertically attaching a metal tube of 2.4-cm diameter to the flow divider, pouring coarse sand into the metal tube, and slightly compacting the sand by knocking on the tube's side-walls. The loam was packed similarly as for the homogeneous columns, leading to a comparable bulk density (Table 1). After both the matrix and PFP domains were established, the metal tube was carefully pulled out of the column. The matrix and PFP domains were thus in direct contact with each other, but without any lateral soil mixing or intrusions between loam and sand, as was verified in horizontal cross sections taken after completion of the experiments. The soil column was slowly saturated from the bottom with deionized water before the first (drainage) experiment.
For infiltration under constant head boundary condition, a tension infiltrometer with infiltration disc of matching diameter (24 cm) was employed. Column outflow through the bottom boundary (consisting of a perforated plate with nylon membrane) was monitored separately for the matrix and the PFP. Pressure heads within the PFP were measured with custom mini-tensiometers (porous ceramic cup: 0.2-cm diam., 1.5-cm length) at depths of 5, 20, 35, 50, and 65 cm. At the same depths, another set of tensiometers with a 0.6-cm-diameter ceramic cup was installed in the soil matrix at distances of 6 and 2 cm from the macropore wall, respectively.
Corresponding to the mini-tensiometers, 8-cm-long TDR probes were used to measure the matrix water contents at the same depths of 5, 20, 35, 50, and 65 cm. Conventional TDR probes are too large to measure water contents in the PFP. Nissen et al. (1998) presented coiled TDR probes with smaller dimensions that allow measurements along short distances. We found these to be more suitable to measure water contents in the PFP, even if they are inherently less accurate than TDR probes. Four TDR coil probes (0.3-cm diam., 2-cm length) were built according to the design of Nissen et al. (1998) and used to monitor the water content at depths of 20 and 35 cm inside the sand-filled PFP and within the loam matrix at a distance of 2 cm from the PFP. The TDR and TDR coil probes were calibrated separately to measure water contents in the matrix loam and PFP sand.
Three preferential flow experiments were selected for this study:
- Infiltration under zero pressure head into the loam column with a medium-sand PFP.
Based on tensiometer readings, the initial hydrostatic pressure head profile increased linearly from –100 cm at the top to –20 cm at the bottom of the soil column. A constant pressure head of –20 cm was imposed at the bottom boundary.
- Infiltration at zero applied pressure head into the loam column with a coarse-sand PFP for somewhat drier initial conditions.
The hydrostatic pressure head profile varied between –150 cm at the top and –70 cm at the bottom. The bottom outlets of the PFP and the matrix were open to the atmosphere.
- Drainage of the initially fully saturated loam column with a coarse-sand PFP.
A pressure head of –10 cm was applied at the lower boundary until 138 min, when it was decreased to –45 cm.
Field Bromide Displacement Tests
For demonstrating the application and suitability of our parameter estimation schemes, we used a field Br– tracer data set from Lennartz et al. (1999). The experimental field site Bokhorst (northern Germany) had a monitored area of 5010 m2. Tile lines were installed at a depth of 1 m and at a spacing of 11 to 14.5 m. The loamy soil at the site was classified as a poorly drained Dystric Gleysol. The 0.5-ha experimental area was identified as heterogeneous in terms of particle size distribution. While the clay content decreased with soil depth from 18% to about 10% at the 1-m depth, bulk density increased from 1.46 to 1.9 g cm–3. The high clay content of the soil surface layer induced cracks and fissures during dry periods. The conventional tilled experimental plot was under a 3-yr crop rotation with winter rape (Brassica napus L. var. napus) (1991/1992), winter wheat (Triticum aestivum L.) (1992/1993), winter barley (Hordeum vulgare L.), (1993/1994), and winter rape (1994/1995). In three Br– tracer tests conducted between 1991 and 1994, preferential flow was found to be a major factor controlling solute transport. Bromide tracer was applied once every year in late fall as a KBr solution, and was monitored in the drainage discharge during the following 5 mo. Applied mass was comparable among the three tracer tests (ranging between 16.75 and 19.3 kg). Readers are referred to Lennartz et al. (1999) for further details about the site, soil, and the tracer experiments.
Numerical Model Simulations
Flow in Homogeneous Soil Columns
The water flow experiments were described by fitting the numerical solution of the Richards' equation as implemented in HYDRUS-1D (imnek et al., 1998) to the observations of inflow; outflow; pressure heads at the 10-, 20-, and 30-cm depths; and water contents at the 10- and 30-cm depths. The resulting van Genuchten (1980) soil hydraulic parameters represented the matrix (loam) and two PFP materials (medium sand and coarse sand). The residual water content, 
r (L3 L–3), was set to zero for the medium and coarse sand. This assumption did not adversely affect the description of water flow in the relatively wet water content range as operated in the water flow experiments. To facilitate the inverse parameter identification for loam and medium sand, the saturated water content, s (L3 L–3), was fixed at the measured value.
Preferential Flow in a Soil Column
Dual-permeability models describe preferential flow using one equation for flow in the matrix and one for flow in the fracture pore system or PFPs, with the two flow regions coupled by a water transfer term. The DPM of Gerke and van Genuchten (1993a) as implemented into the modified HYDRUS-1D (imnek et al., 2003) uses the following two mixed-type Richards equations to describe transient water flow in fractures (Eq. [1a]) and the matrix (Eq. [1b]):
 | [1a] |
 | [1b] |
where the subscript "f" defines a property of the fracture pore system or PFP, the subscript "m" represents the matrix, is the water content (L3 L–3), h is the pressure head (L), K denotes the hydraulic conductivity function (L T–1), wf is the fraction of total soil occupied by the fracture system (0 < wf < 1), and 
w is the water transfer term (T–1)
 | [2] |
where ß is a dimensionless geometry coefficient, a is the characteristic length of the matrix structure (L), Ka is the hydraulic conductivity of the fracture–matrix interface (L T–1), and 
w is a dimensionless scaling coefficient for which a value of 0.4 was suggested (Gerke and van Genuchten, 1993b). Because of the absence of depositional coatings at the PFP–matrix interface in our repacked soil columns, we assumed that Ka in the transfer term can be described as
 | [3] |
The geometry factor ß was calculated assuming hollow cylindrical PFP geometry according to Gerke and van Genuchten (1996) as
 | [4] |
where the parameter 
was calculated as
 | [5] |
in which b is the cylindrical PFP radius. The volume fraction of the PFP, wf, was calculated as
 | [6] |
The numerical grid consisted of 201 elements with a vertical nodal spacing of 0.4 cm to represent the 80-cm-long column. The values for mass-transfer related DPM parameters were derived from the experimental PFP–matrix geometry (a = 11-cm matrix mantle radius, b = 1.2-cm PFP radius, = 10.167, ß = 1.0685, and wf = 0.009675).
The initial condition was defined by the hydrostatic pressure head profile interpolated between initial tensiometer readings at various depths. The lower boundary in the model was represented by a seepage face boundary condition, with no flow when the boundary was unsaturated, and a zero (atmospheric) pressure head once saturation was reached. For those experiments where suction was applied to the lower boundary, the seepage face boundary condition was modified by using the negative value of applied suction as the threshold value between zero flux and constant pressure head.
Field-Scale Preferential Flow and Bromide Transport
For simplicity, the entire 0.5-ha tile-drained field was represented by a one-dimensional vertical single-layered model domain, since assuming several layers or accounting for two-dimensional heterogeneity was impractical for inverse DPM simulation because of the increasing number of parameters. Soil hydraulic and solute transport DPM parameters of the one-dimensional homogeneous domain were assumed to represent effective soil properties lumped over heterogeneities in soil texture, structure, and effects of tile-induced flow-field.
Preferential water flow was described with Eq. [1] through [3]. However, since the parameters a, Ka, ß, and w in Eq. [2] are difficult to define on the field scale, they were lumped into a single water transfer coefficient, 
w (T–1 L–1), which yields the following simplified expression for the water transfer term:
 | [7] |
Convection–dispersion equations simulate transport in both the fracture and matrix regions, again using the DPM of Gerke and van Genuchten (1993a) as implemented into the modified HYDRUS-1D (imnek et al., 2003), as follows:
 | [8a] |
 | [8b] |
where for the fracture and matrix regions, c is the solute concentration (L3 M–1), q is the volumetric flux density (L T–1), and D is the dispersion coefficient (L2 T–1) described with the following equation:
 | [9a] |
 | [9b] |
where once again for the fracture and matrix regions, D0 represents the molecular diffusion (L2 T–1), 
signifies the dispersivity (L), and s denotes the saturated water content. The rate of solute mass transfer between the fracture and matrix regions, s, is given as the sum of diffusive and convective fluxes:
 | [10] |
where s is a first-order diffusive solute mass transfer coefficient (T–1), and c* is equal to cf for w > 0 and cm for w < 0.
The initial condition for water flow was assumed as a hydrostatic profile between zero pressure head at drainage depth (100 cm) and –100 cm at the top. Measured rainfall rates less the actual evapotranspiration were applied as the upper boundary condition. A seepage face boundary condition was prescribed for water flow at the lower boundary. The initial Br– concentration was assumed to be zero in 1991 and 1993. In 1994, resulting from the prior tracer tests, an initial concentration of 4 mg L–1 was assumed throughout the soil. The top boundary condition for solute transport was represented by the prescribed Br– concentration of the applied KBr solution, while a zero Br– concentration gradient was prescribed at the lower boundary.
Inverse Parameter Identification
Preferential Flow in a Soil Column
For describing preferential water flow in the soil column using Eq. [1] to [6], the DPM requires 13 parameters (i.e., m,r, m,s, 
m, nm, Km,s, f,r, f,s, f, nf, Kf,s, a, wf, ß), of which three parameters related to water transfer between domains are known (a, wf, ß). To additionally reduce the number of unknowns, the residual water content for the fracture pore system, f,r, was fixed at zero. This limited the total number of unknown water flow parameters to nine. The DPM model parameters were estimated by minimizing an objective function, 
, (imnek et al., 1998) representing the deviations between measured and calculated variables as follows:
 | [11] |
where m is the number of different sets of measurements, n represents the number of observations in a particular measurement set, Oj(z,ti) represents observations at time ti for the jth measurement set at location (z), Ej(x, ti, b) are the corresponding estimated space–time variables for the vector b of nine optimized van Genuchten (1980) parameters (m,r, m,s, m, nm, Km,s, f,s, f, nf, Kf,s) and wi,j are weighting factors that here were manually calibrated in repeated inverse simulation runs until similar contributions of all observation sets to the weighted least squares (Eq. [11]) of residuals between observations and model estimates (Eq. [11]) were achieved. This ad hoc procedure gave more evenly distributed deviations between individual data sets and model results than using other weighting options for wi,j, such as the ratio of means or standard deviations of data sets. However, improved automated weighting procedures should be developed in future. The Levenberg-Marquardt algorithm was applied to minimize the objective function (Eq. [11]) (imnek et al., 1998).
For the case of drainage of the soil column initially completely saturated up to the top, a modification of the above procedure became necessary to properly describe the gradual, noninstantaneous decrease of positive initial pressure heads (see results). At the bottom of the matrix and PFP domains, a thin layer corresponding to the bottom perforated-plate thickness (L = 0.4-cm height) with a constant hydraulic resistance factor of R = L/Ks was assumed. The factor R accounted for any potential hydraulic resistance of the bottom plate with nylon membrane during the early phase of drainage, characterized by positive pressure heads and large hydraulic gradients across the lower boundary. For technical reasons, R could not be measured below matrix and PFP, and thus for the drainage case, R for matrix and PFP had to be inversely estimated. Allowing for different magnitudes of R at the bottom of matrix and PFP was justified by the fact that the perforated-plate had a higher perforation density below the PFP than below the matrix.
Three different approaches were compared for estimating the parameters of the DPM. First, the inverse-lumped approach relied on observations of pressure heads in the matrix at 5-, 35-, and 65-cm depths, water contents in the matrix at 35-cm depth, and cumulative infiltration (across top boundary) or drainage (across bottom boundary) fluxes. Second, the inverse-local approach additionally utilized observations of region-specific (matrix and PFP) outflow, that is, data that were not available in earlier column studies or field experiments. Finally, in the forward approach we studied if independently determined soil hydraulic parameters for loam, medium sand, and coarse sand could be successfully used for forward DPM simulation of the preferential flow data in the repacked soil columns with various matrix–PFP compositions.
Field-Scale Preferential Flow and Bromide Transport
Two inverse DPM approaches were compared to simulate the field experiments of Lennartz et al. (1999). In the sequential DPM approach, eight hydraulic parameters (i.e., m, nm, Km,s, f,s, f, nf, Kf,s, w) were estimated based on observations of tile-drainage outflow (both in terms of rate and cumulative values), followed by the estimation of one to three solute transport parameters for different scenarios (s; or s, f; or s, f, m) based on the observed Br– concentrations in tile-drainage effluent. In the simultaneous DPM approach, all eight hydraulic and one to three solute parameters were estimated at the same time using both tile-drainage and Br– concentrations. In both approaches, the parameters were estimated by minimizing the objective function (Eq. [11]) as described above for simulations of soil column tests. The following parameters were fixed: the volumetric portion of the fracture pore system (wf = 0.05), as estimated from dye-stained soil fractions (Köhne, 1999), the residual water content of the fracture pore system (f,r = 0), residual and saturated matrix water contents (m,r = 0.078, m,s = 0.43) of a loam as predicted by the neural network prediction option in HYDRUS-1D (imnek et al., 1998), and the matrix dispersivity at an average magnitude (m = 0.5 cm) found in Br– displacement tests conducted with Bokhorst soil columns (Meyer-Windel, 1998).
Goodness-of-Fit
The goodness-of-fit of the three parameter modeling schemes for preferential flow (inverse-local, inverse-lumped, and forward) was evaluated using the modified index of agreement (MIA) as given in Legates and McCabe (1999). For each measurement set j, the index MIAj was calculated as follows:
 | [12] |
The MIA (here: MIAj) varies between 0.0 to 1.0, with higher values indicating better agreement between the model and observations, similar to the interpretation of the coefficient of determination, R2. Compared with R2, MIA is a more reliable statistical measure since it is more sensitive to differences in the observed and model simulated means and variances and is less sensitive to extreme values (Legates and McCabe, 1999). There is no statistical basis to decide which MIA value is a threshold separating "inaccurate" from "accurate" simulations. Here we assumed all column-scale water flow simulations with MIA > 0.75 to be "accurate." However, this value would be too rigorous for assessing field-scale preferential solute transport simulations, for which other criteria are usually applied, such as the percentage of the simulation results within a certain deviation from observations (e.g., Larsson and Jarvis, 1999a, 1999b). Hence, a common MIA threshold value to separate accurate from inaccurate hydraulic and transport simulations at field scale is not applicable. However, visual comparison of the agreement between model result and data may provide additional information (Green and Stephenson, 1986). Therefore, for the field-scale simulations, MIA values were compared only by relative ranking, and additionally two subjective binary (yes/no) goodness-of-fit measures were employed based on visual comparison of observed and simulated Br– breakthrough curves: (i) match of first preferential Br– peak and (ii) match of breakthrough curve pattern (e.g., peaked vs. smooth, skewed vs. symmetric).
|
RESULTS AND DISCUSSION
|
|---|
Homogeneous Soil Columns
The experimental water flow observations and corresponding inverse simulation results for the Richards' equation using inverse-HYDRUS-1D are documented briefly (Fig. 1
and 2)
. The corresponding estimated van Genuchten's (1980) soil hydraulic parameters that were subsequently used as region-specific parameters in the DPM forward approach are given in Table 2. Figure 1 shows results of infiltration and drainage experiments conducted in the loam column, such as cumulative inflow and outflow (Fig. 1a and 1b), pressure heads (Fig. 1c and 1d) and water contents (Fig. 1e and 1f). Figure 2 presents experimental and inverse simulation results for infiltration into the medium-sand column and for drainage of the initially saturated coarse-sand column at three applied bottom suctions of 0, –30, and –60 cm in sequence. For technical reasons, water content measurements could not be retrieved during the drainage experiment in the coarse sand except for the saturated initial water content (Fig. 2f). All simulations were accurate according to the MIA values exceeding 0.75.
View larger version (32K):
[in this window]
[in a new window]
|
Fig. 1. Experimental data and inverse HYDRUS-1D simulation results for the loam column, infiltration (left) and drainage (right): (a, b) cumulative infiltration and outflow, (c, d) pressure heads at the 10-, 20-, and 30-cm depths, and (e, f) water contents at the 10- and 30-cm depths.
|
|
View larger version (33K):
[in this window]
[in a new window]
|
Fig. 2. Experimental data and inverse HYDRUS-1D simulation results for the infiltration into the medium-sand (left) and coarse-sand (right) columns: (a, b) cumulative infiltration and outflow, (c, d) pressure heads at 10-, 20-, and 30-cm depths, (e, f) water contents at the 10- and 30-cm depths.
|
|
Infiltration into the Loam Column with Medium-Sand Preferential Flow Path
The case of infiltration into the loam column containing a PFP of medium sand was selected to represent an example of relatively small hydraulic nonequilibrium or low degree of preferential flow. Figure 3
shows experimental and simulated results for the infiltration into the loam column with medium-sand PFP. The observed total infiltration of 1067 cm3 of water (Fig. 3a) exceeded the total outflow of 665 cm3 (Fig. 3b) by 400 cm3; that is, there was an increase in water storage of about 40%. The cumulative region-specific outflow via the matrix of 380 cm3 (Fig. 3c) and via the medium-sand PFP of 230 cm3 (Fig. 3d) was equivalent to a 60-fold local water flux out of the PFP as compared with the matrix, considering the 100:1 volume ratio of matrix to PFP. The corresponding larger Ks ratio of PFP to matrix materials of 350 (609/1.7, see Table 2) suggests that outflow out of the PFP was considerably reduced by lateral water transfer from the PFP into the matrix.
View larger version (23K):
[in this window]
[in a new window]
|
Fig. 3. Infiltration into the loam column with a preferential flow path (PFP) made of medium sand. Comparison of experimental cumulative flow data and dual-permeability model simulations using independently estimated parameters (forward), or inversely identified parameters based on data of pressure heads and water contents in the matrix, and domain-specific (inverse-local) or bulk-soil related (inverse-lumped) outflow. (a) Infiltration, (b) outflow, (c) outflow from the matrix, and (d) outflow from the PFP.
|
|
Both the inverse-local (Table 3) and inverse-lumped (Table 4) DPM approaches accurately matched the bulk soil related observations of cumulative infiltration, qin (Fig. 3a), and cumulative outflow, qout (Fig. 3b), as revealed by the large MIA values (0.95–0.97, Table 5). The forward DPM simulation also accurately (MIA = 0.92) predicted the measured outflow (Fig. 3b) without any calibration, while the MIA for infiltration (0.72) was only just below the assumed threshold value of 0.75 for an accurate simulation. The region-specific outflow observations for the matrix (Fig. 3c) and PFP (Fig. 3d) regions were accurately described by inverse-local (MIA 
0.95) and forward DPM (MIA > 0.75) simulations, while the inverse-lumped approach overestimated outflow out of the matrix, qm,out, and compensated this overestimation by a corresponding underestimation of outflow out of the PFP, qf,out (see Fig. 3c and 3d), yielding MIA values for qm,out and qf,out far below 0.75 (Table 5).
View this table:
[in this window]
[in a new window]
|
Table 5. Modified index of agreement, MIA (Eq. [12]), for differences between various model results and experimental observations.
|
|
Figure 4
compares experimental and simulation results of pressure heads in the matrix and in the PFP during the first 1.5 h at the 5- (Fig. 4a), 35- (Fig. 4b), and 65-cm (Fig. 4b) depths. The pressure heads in the matrix and the PFP were almost similar for a particular depth, revealing the low degree of hydraulic nonequilibrium. Both the inverse-local and the forward DPM gave a succession of steeper pressure head increases with depth, which suggested a somewhat sharper infiltration front traveling through the column than would be assumed for more gradually increasing pressure head data. Furthermore, the forward simulation gave pressure heads in the PFP that increased ahead of those in the matrix. In the inverse-local simulation, this was the case at 5 cm below the top of the column, while at the 65-cm depth the trend was reverse (i.e., matrix pressure heads increased earlier than PFP). The reason for this discrepancy is the difference between region-specific hydraulic conductivity functions: for inverse-local, the crossover-point where the PFP hydraulic conductivity falls below that of the matrix was at h = –19 cm, while for forward it was at h = –120 cm (Table 2). Consequently, in the inverse-local simulation the matrix had a higher hydraulic conductivity in part of the unsaturated pressure head range observed in the experiment (–40 cm < h < 0).
View larger version (28K):
[in this window]
[in a new window]
|
Fig. 4. Same as Fig. 3, but showing pressure heads in matrix and preferential flow path (PFP) at depths of (a) 5, (b) 35, and (c) 65 cm.
|
|
Although the inverse-local and the forward DPM simulations were in qualitative agreement with pressure heads (Fig. 4), the MIA values below 0.75 indicated that pressure heads were not accurately described by any approach. One reason for the deviations between model simulations and data may be the inherent inaccuracy of the first-order water exchange term (Köhne et al., 2004). Additionally, region-specific pressure head measurements may not have been as accurate as flow measurements. However, we note that matrix pressure heads (and matrix water contents shown below) in this study are regarded as auxiliary information only, the primary information being the infiltration and outflow, since outflow was also the only hydraulic information available and used in the field-scale simulations.
Water contents in the matrix and in the PFP at the 35-cm depth are shown in Fig. 5
. The matrix water contents stayed close to saturation without showing any clear reaction to infiltration. The water contents in the PFP increased even more gradually than pressure heads at the 35-cm depth (see Fig. 4 and 5). However, it cannot be entirely ruled out that part of the TDR coil probe was not placed in the PFP, which would have resulted in a measured water content averaged over the PFP and adjacent matrix. The forward simulation resulted in the best prediction of water contents in the PFP. The inverse-local approach produced a relatively sharp PFP water content rise that was not characteristic of the observed PFP water contents (Fig. 5). The assumption of zero residual water content of the medium-sand PFP, made to reduce the number of unknowns, may not have been realistic. The small MIA values (Table 5) were in sharp contrast with the inverse-local approach clearly matching the matrix water contents (Fig. 5), such that the inadequacy of the MIA for this type of data scattering around an almost constant level became evident.
The inverse-lumped DPM simulation results (not shown) generally showed slightly stronger deviations from the pressure head and water content data than the inverse-local and the forward DPM approaches. The results were omitted in Fig. 4 and 5 since the inverse-lumped approach already failed to describe the region-specific outflow.
Infiltration into the Loam Column with Coarse-Sand Preferential Flow Path
The infiltration into the loam column containing a PFP filled with coarse sand constitutes an example of relatively strong hydraulic nonequilibrium or pronounced preferential flow. Figure 6
shows the observed and simulated cumulative flow at the top and bottom of the column. Infiltration (Fig. 6a) and outflow out of bulk soil (Fig. 6b) as well as outflow out of the PFP (Fig. 6d) were much faster, while outflow out of the matrix (Fig. 6c) was less than for the column with medium-sand PFP (Fig. 3). In this way, as expected, fast preferential flow in the coarse-sand PFP was relatively less impeded by transfer into the matrix than flow in the medium-sand PFP.
View larger version (24K):
[in this window]
[in a new window]
|
Fig. 6. Infiltration into loam column with preferential flow path (PFP) made of coarse sand. Comparison of experimental cumulative flow data and dual-permeability model simulations using independently estimated parameters (forward), or inversely identified parameters based on data of pressure heads and water contents in the matrix, and domain-specific (inverse-local) or bulk-soil related (inverse-lumped) outflow. (a) Infiltration, (b) outflow, (c) outflow from the matrix, and (d) outflow from the PFP.
|
|
As evident from MIA values (Table 5) and visual assessment (Fig. 6), the inverse-local DPM simulations once again gave an accurate description of both global (bulk soil related) and local (PFP) flows, and the inverse-lumped DPM approach provided a similar match of global flow data, while PFP flow data were matched less. For the cumulative matrix outflow exhibiting a very flat increase, the low MIA value of 0.07 suggested failure even for the inverse-local simulation, which according to visual evaluation was only slightly overestimating the outflow (Fig. 6). Again, as already found above for the water content data, a conclusion about goodness-of-fit of one subset of observations should not be drawn from a statistical measure like the MIA alone. We note also that R2 values (not shown) were consistently higher than MIA values, which shows that the MIA is a more rigorous statistical measure. For example, an R2 value of 0.74 corresponds to the above-mentioned MIA value of 0.07.
We did not include results of the forward DPM simulation because independently obtained parameters were affected by finger flow phenomenon in the homogeneous coarse-sand column, and we inferred that these parameters would not be representative of the PFP.
The observed and simulated pressure heads in matrix and coarse-sand PFP are shown in Fig. 7
for the depths of 5 (Fig. 7a) and 35 cm (Fig. 7b). The pressure heads showed more pronounced differences between the (coarse-sand) PFP and matrix than was observed between the matrix and a PFP made of medium sand (Fig. 4). This finding of stronger pressure head nonequilibrium can be explained with the larger contrast between hydraulic properties of loam and coarse-sand PFP than between hydraulic properties of loam and medium-sand PFP. Although within the pressure head range shown in Fig. 7 the inverse-local approach apparently better described matrix pressure head observations than the inverse-lumped approach, both gave accurate representations according to similar MIAs of about 0.78 (Table 5). Moreover, both inverse-local and inverse-lumped simulations slightly overestimated water contents in the matrix at the 5-cm depth and underestimated them at the 35-cm depth (Fig. 8
).
Drainage of the Loam Column with Coarse-Sand Preferential Flow Path
Finally, the drainage of the initially entirely saturated loam column containing a PFP made of coarse sand was conducted to see if findings for the infiltration cases could be generalized to drainage flow conditions as well. Findings are presented for a two-step drainage outflow experiment. Figure 9
reveals that cumulative water outflow of the bulk soil (Fig. 9a), soil matrix (Fig. 9b), and PFP (Fig. 9c) increased faster initially and showed slow response to the step change of suction at the lower boundary from 10 to 45 cm at 138 min. The inverse-local DPM approach again provided an accurate match of both lumped outflow and local outflow components, yielding MIA values >0.75 (Table 5). The inverse-lumped DPM matched lumped outflow even better according to the MIA (0.97 vs. 0.92), but once again failed to describe local outflow components of the matrix and PFP (MIA < 0.75, Table 5). The forward DPM approach underestimated outflow out of the matrix (Fig. 9b) and grossly underestimated outflow out of the coarse-sand PFP (Fig. 9c). As a result of this behavior, outflow out of bulk soil was also underestimated by the forward DPM approach (Fig. 9a), and low MIA values were obtained (Table 5).
View larger version (26K):
[in this window]
[in a new window]
|
Fig. 9. Drainage from the loam column with the preferential flow path (PFP) made of coarse sand. Comparison of experimental cumulative outflow data and dual-permeability model simulations using independently estimated parameters (forward), or inversely identified parameters based on data of pressure heads and water contents in the matrix and domain-specific (inverse-local) or bulk-soil related (inverse-lumped) outflow. (a) Outflow, (b) outflow from the matrix, and (c) outflow from the PFP.
|
|
The pressure heads in the matrix and in the PFP at the 5-cm depth (Fig. 10a
) decreased starting from slightly negative values. Pressure heads at the 65-cm depth started from initially positive values, decreased within 6 min (PFP) and 20 min (matrix) to negative values, and showed a reaction to the suction step increase at 138 min (Fig. 10b). Differences between pressure heads in PFP and matrix regions were less than for the infiltration case (compare Fig. 7 and 10), which reveals that the degree of pressure head nonequilibrium during preferential flow not only depends on the difference of hydraulic properties between PFP and matrix, but also on the flow direction (infiltration vs. drainage). The inverse-local and inverse-lumped simulations described pressure heads accurately (MIA above 0.75; see Fig. 10) when accounting for a hydraulic resistance at the bottom. Pressure heads simulated with the forward DPM approach did not match with the observations (Fig. 10, Table 5).
View larger version (26K):
[in this window]
[in a new window]
|
Fig. 10. Same as Fig. 9, but showing pressure heads in the matrix and in the preferential flow path (PFP) at depths of (a) 5 and (b) 65 cm.
|
|
Water contents at the 35-cm depth showed hardly any reaction with time for the matrix, but a strong decrease in the water content of the PFP (Fig. 11
) occurred during the onset of drainage. According to visual evaluation, the inverse models gave a satisfactory representation of the PFP water contents and excellent description of matrix water contents not reflected in the MIA (Table 5), while the forward simulation results deviated relatively more from the experimental observations (Fig. 11).
Dual-Permeability Model Parameters and Their Uncertainty
Model parameters and their 95% confidence intervals for the forward, inverse-local, and inverse-lumped approaches are listed in Tables 2, 3, and 4. Results for steady-state Ks and static retention curve measurements are given for comparison (Table 2). The differences between parameter magnitudes for the various approaches reflect the observed deviations among simulation results, which will be discussed for the Ks values as an example. For the infiltration experiments, Ks for loam was lower for the inverse-local approach, with values of 3.9 and 3.1 cm d–1 for columns with coarse- and medium-sand PFPs, respectively (Table 3), than corresponding loam Ks values for the inverse-lumped approach of 8.4 and 11.6 cm d–1, respectively (Table 4). The large matrix Ks values obtained with the inverse-lumped approach resulted in the overestimation of the observed cumulative matrix outflow (Fig. 3 and 6). Accordingly, the values obtained in the inverse-local approach compared better, with Ks of 1.7 cm d–1 obtained for infiltration into the homogeneous loam column. Assuming Ks equal to 1.7 cm d–1 in the forward DPM approach still produced accurate prediction of matrix outflow (Fig. 3). On the other hand, once again for the infiltration experiments, the inverse-local approach gave higher Ks estimates for mediums and coarse-sand PFPs of 550 and 4343 cm d–1, respectively (Table 3), than the inverse-lumped approach (138 and 3380 cm d–1, Table 4), which for the inverse-lumped approach resulted in the underestimation of the observed cumulative PFP outflow (Fig. 3 and 6). Consistent with these observations, for the medium-sand PFP the Ks of 550 cm d–1 of the inverse-local approach was closer to Ks obtained for infiltration in the medium-sand column of 609 cm d–1, than Ks of 138 cm d–1 for the inverse-lumped approach. Moreover, Ks confidence limits obtained for infiltrations in matrix and PFP were consistently smaller for the inverse-local than inverse-lumped approach (Tables 3, 4). Furthermore, the correlation, r, between Ks,f and Ks,m values in PFP and matrix was small for the inverse-local approach (r = –0.33), but an almost perfect negative correlation (r = –0.99) was obtained for the inverse-lumped approach (Table 6). This highly negative correlation of the inverse-lumped approach means that a higher Ks,f estimate is compensated by a lower Ks,m value and vice versa. In other words, Ks,f and Ks,m cannot be inversely identified at the same time from bulk soil water flow observations. In total, four "large" correlations (r > 0.75) were identified for the inverse-lumped approach, but only two values of r > 0.75 for the inverse-local approach, between and n of both matrix and PFP (Table 6). However, in the case of empirical curve-fit parameters without much physical meaning, such as and n, correlations may be acceptable if only the shape of the estimated hydraulic functions is of interest, which may be similar for different combinations of such correlated parameters.
View this table:
[in this window]
[in a new window]
|
Table 6. Correlation, r, between parameters for inverse dual-permeability model simulations (inverse-lumped and inverse-local approaches) of infiltration into loam with coarse-sand PFP.
|
|
All the above comparisons together suggest that the inverse-local approach not only accurately described preferential water flow, but that the resulting inverse Ks estimates (and other parameters except and n) were also generally well identified and physically meaningful, while all this was not the case for the inverse-lumped approach. This is a key finding of this study, the practical implications of which are shown below for DPM prediction of field-scale preferential solute transport based on tile-drainage hydrographs.
For the drainage experiments and both inverse DPM approaches, Ks,f and Ks,m and s,m were larger than for infiltration (Tables 3 and 4). This discrepancy can be explained when considering the two wetting fronts in the matrix during preferential infiltration, that is, one from the top and one from the side by water transfer from the PFP. Hence, beneath the vertical wetting front in the matrix, there were coexisting wet (adjacent to PFP) and dry (far from PFP) regions, which resulted in lower Ks and s estimates, as compared with values for water drainage coming from the total matrix cross section. Additionally, lateral water transfer during infiltration may have entrapped air below the matrix wetting front, also resulting in lower Ks and s estimates for the matrix. In summary, the above findings are a manifestation of enhanced soil hydraulic hysteresis during preferential flow. Indications of hysteresis during nonequilibrium flow have been found before (e.g., imnek et al., 2001; Köhne and Mohanty, 2005). The relations between Ks values found for the various approaches were not as clear for drainage as was found for infiltration, which we attribute to the difficulties in modeling the observed gradual positive pressure head dissipation, requiring the estimation of lower boundary resistivity factors for matrix and PFP. Moreover, in the drainage case, the hydraulic resistance factor greatly increased the uncertainty associated with the Ks values of matrix and PFP, as apparent from large confidence intervals (Tables 3 and 4).
For comparison, the magnitudes of Ks measured at steady state for loam and medium and coarse sand were higher than measured under transient flow conditions in homogeneous columns, but the corresponding 95% Ks confidence limits were still overlapping. The s values obtained in static measurements were generally between those obtained for infiltration and drainage experiments (Tables 2–
4).
Those parameter values that were strongly correlated are no longer amenable to interpretation and comparison, while the shapes of hydraulic functions, which may be similar for different values of correlated parameters, can still be compared across model approaches. Figure 12
compares water retention and hydraulic conductivity functions obtained for homogeneous columns and steady-state or static conditions with those obtained for the matrix and PFP. For medium sand, inverse-local and-lumped produced similar PFP retention functions (Fig. 12a) and hydraulic conductivity functions (Fig. 12c), but these functions were higher close to saturation and lower in the drier range than those obtained for infiltration into the medium-sand column. For coarse sand, the inverse-lumped water retention (Fig. 12b) and hydraulic conductivity (Fig. 12d) functions were vastly different, which was unrealistic considering that retention functions obtained by static and infiltration measurements were quite similar (Fig. 12b). For the loam, retention functions of all approaches were almost parallel and revealed that estimated different r values (and maybe n and values) obtained in different approaches are not necessarily meaningful and probably cannot be used outside the experimental pressure head range near saturation. The different (convex or concave) shapes of the loam hydraulic conductivity function for infiltration and drainage were better reproduced by the inverse-local approach, disregarding differences in the absolute values (Fig. 12f).
View larger version (40K):
[in this window]
[in a new window]
|
Fig. 12. Soil hydraulic van Genuchten functions fitted to the soil water retention curve measured at equilibrium (static) and obtained by fitting the numerical solution of Richards' equation to water flow measurements performed on homogeneous soil materials (infiltration and drainage), and obtained by fitting the dual-permeability model to flow measurements conducted with heterogeneous soil systems consisting of loam (matrix) and medium or coarse sand (preferential flow path, PFP): medium sand (a) water content and (b) hydraulic conductivity, coarse sand (c) water content and (d) hydraulic conductivity, and loam (matrix) (e) water content and (f) hydraulic conductivity.
|
|
Simulation of Field Tracer Tests
Results for simulated drain discharge rates and cumulative drain discharge using the sequential and simultaneous DPM simulation schemes are compared in Fig. 13
with corresponding experimental observations published by Lennartz et al. (1999). For completeness, the natural daily rainfall used as the upper boundary condition is provided as well (Fig. 13). The drain discharge rates closely reflect the rainfall pattern and were reproduced by both inverse DPM schemes, the simultaneous scheme yielding somewhat higher and narrower discharge peaks than the sequential scheme. The resulting MIA values were only slightly larger for the sequential approach (Table 7). Cumulative drain discharge is reproduced equally well with both schemes, according to MIA values around 0.97, except for 1994/1995 where the sequential approach yielded an MIA of 0.98, as compared with 0.91 for the simultaneous approach.
View larger version (30K):
[in this window]
[in a new window]
|
Fig. 13. Rainfall, drain discharge flux, and cumulative drain outflow as observed at the Bokhorst field site for three runoff seasons between 1991 and 1995 (data: Lennartz et al., 1999), and inverse dual-permeability model (DPM) simulation results for the sequential and simultaneous approaches.
|
|
View this table:
[in this window]
[in a new window]
|
Table 7. Statistical goodness-of-fit: modified index of agreement, MIA (Eq. [12]), for differences between model results and observations (Lennartz et al., 1999) of drain discharge flux (qout), cumulative drain discharge (qcumout), and Br– concentrations in tile drainage water [c(Br–)]; and success (Yes) or failure (No) of Br– transport simulation based on visual comparison of observed and simulated Br– breakthrough curves (Fig. 14 and 15); match of first preferential Br– peak (First peak); and match of general breakthrough curve pattern with regard to shape and secondary peaks (Pattern).
|
|
View larger version (34K):
[in this window]
[in a new window]
|
Fig. 14. Bromide tracer concentrations as observed in drain outflow at the Bokhorst field site for three runoff seasons between 1992 and 1995 (data Lennartz et al., 1999) and inverse dual-permeability model (DPM) simulation results for the sequential and simultaneous approaches. The only parameter fitted in the solute transport equation was the first-order solute transfer coefficient.
|
|
View larger version (28K):
[in this window]
[in a new window]
|
Fig. 15. Same as Fig. 14, but in the solute transport equation, both the first-order transfer coefficient and the fracture dispersivity (1994/1995: additionally the matrix dispersivity) were fitted.
|
|
The observed and simulated Br– concentrations are shown in Fig. 14
for the case with dispersivities in matrix, m, and PFPs, f, fixed at 0.5 cm as derived from Br– displacement experiments using soil columns (Meyer-Windel, 1998). The estimated soil hydraulic and transport DPM parameters are listed in Table 8. In all 3 yr, there was a recurrent pattern of rapid Br– concentration increase to the main preferential Br– peak directly following the first rainfall event after solute application, and secondary smaller Br– concentration peaks in response to later rainfall events. Only toward the end of the 1994/1995 season, was there a slight increase in Br– concentration considered the matrix Br– peak (Fig. 14). Strikingly, only the simultaneous DPM approach consistently reproduc
TOP
ABSTRACT
INTRODUCTION
