Published online 23 August 2007
Published in Vadose Zone J 6:651-667 (2007)
DOI: 10.2136/vzj2007.0033
© 2007 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
ORIGINAL RESEARCH
Two-Dimensional Dual-Permeability Analyses of a Bromide Tracer Experiment on a Tile-Drained Field
Horst H. Gerkea,*,
Jaromir Dusekb,
Tomas Vogelb and
J. Maximilian Köhnec
a Institute for Soil Landscape Research, Leibniz-Centre for Agricultural Landscape Research (ZALF), Eberswalder Strasse 84, D-15374 Müncheberg, Germany
b Dep. of Hydraulics and Hydrology, Faculty of Civil Engineering, Czech Technical Univ., Prague, Czech Republic
c Faculty for Agricultural and Environmental Sciences, Soil Physics and Resources Protection, Univ. of Rostock, Justus-von-Liebig-Weg 6, D-18059 Rostock, Germany
* Corresponding author (hgerke{at}zalf.de).
All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.
Received 10 February 2007.
 |
ABSTRACT
|
|---|
Preferential flow has been hypothesized as an important factor for chemical leaching from tile-drained agricultural fields with structured soils originating from glacial till sediments. Previous studies showed that one-dimensional single-porosity models (1D-SPM) failed and that one-dimensional dual-permeability models (1D-DPERM) were limited in explaining both Br leaching and residual Br distribution, although tile water outflow peaks could somehow be reproduced. The objective of this paper was to analyze the tile outflow and leaching patterns using a two-dimensional (2D)-DPERM and a standard 2D-SPM for comparison. Flow and transport were simulated in a 2D vertical cross-section of 5.9 m length and 2 m depth using previously tested parameters. Simulated drainage rates and Br-effluent concentrations were made comparable with collector data from a field experiment by weighing results for irrigated and nonirrigated plots according to their area fractions. The 2D-DPERM simulations for surface application of Br in dissolved form in both domains overestimated the observed initial outflow concentration peaks, in contrast to closer approximation of observations assuming Br application in the soil matrix domain only. The simulated 2D mass transfer rate distribution showed most intensive exchange between domains near the water table and in the topsoil. Results from the 2D-DPERM analyses suggest that conditions at the soil surface, near the water table, and of the field-scale mixing are significantly affecting leaching patterns, in addition to local nonequilibrium effects. Here, the description of preferential flow toward tile drain could be strongly improved with the 2D-DPERM compared with the 2D-SPM. Further improvements remain challenging with respect to DPERM numerical modeling and field experimentation, with special attention toward soil structure and soil surface conditions.
Abbreviations: 1D, one-dimensional 2D, two-dimensional 3D, three-dimensional BC, boundary condition Br, bromide CDE, convection–dispersion equation DPERM, dual-permeability model DPM, dual-porosity model ET, evapotranspiration MIM, mobile–immobile model PF, preferential flow SM, soil matrix SPM, single-porosity model
 |
INTRODUCTION
|
|---|
Tile-drain systems are widely used for controlling soil moisture of arable fields with temporary high ground or perched water tables to improve conditions for crop growth and soil cultivation. Subsurface drainage not only affects the field hydrology but also bears a potential for chemical leaching together with the tile effluent (Vereecken, 2005; Schelde et al., 2006). Leaching may be enhanced especially by preferential flow in structured soils (Bosch et al., 2000), thereby increasing the risk of pollution of nearby surface water by agricultural chemicals (Stamm et al., 1998; Kohler et al., 2001; Stone and Wilson, 2006). Numerous tracer tests have been performed on tile-drained arable fields to study preferential field-scale transport of solutes (Czapar et al., 1994; Bronswijk et al., 1995; Kung et al., 2000a, 2000b; Jaynes et al., 2001; Stamm et al., 2002) and pesticides (Czapar et al., 1994; Elliott et al., 2000; Kladivko et al., 2001; Zehe and Flühler, 2001; Kronvang et al., 2004). Such experimental studies use the fact that tile drainage and leaching rates integrate all effects of soil spatial variability and soil structure occurring at the field site, which may be viewed as a field-scale lysimeter (Richard and Steenhuis, 1988).
Preferential flow may be identified by the event-driven simultaneous early appearance of tracer and reactive chemicals in tile outflow (e.g., Lennartz et al., 1999; Zehe and Flühler, 2001; Fortin et al., 2002). For a tile-drained experimental field site with structured soils developed in glacial till sediments (Bokhorst, northern Germany), preferential flow showed characteristic patterns of highest tracer concentrations during the first leaching event after application and concentration peaks within the initial phase of subsequent events (Lennartz et al., 1999; Köhne, 1999). Similar concentration patterns in tile effluent were observed in other field tracer experiments in postglacial landscapes (Villholth et al., 1998; Mohanty et al., 1998; Jaynes et al., 2000; Haws et al., 2004, 2006). Attempts to explain observed tile effluent patterns include 1D and 2D numerical models with a single pore domain (single-porosity model [SPM]) and two pore domains (dual-porosity model [DPM]) (Haws et al., 2005; Boivin et al., 2006; Köhne et al., 2006) and transfer function models (Arabi et al., 2006; Gaur et al., 2006), among other approaches (Kohler et al., 2003). The two-domain concept, assuming either one (mobile–immobile model [MIM]) or two flow domains (mobile–mobile or dual-permeability model [DPERM]), has been proposed to describe transient flow and solute transport in structured soil under variably saturated conditions (
im
nek et al., 2003; Gerke, 2006). The DPERM descriptions assume either gravity-driven flow in macropores (Jarvis, 1994; Larsson and Jarvis, 1999; Larsbo et al., 2005) or Darcian flow in a more permeable (macro-) pore system (Gerke and van Genuchten, 1996; Vogel et al., 2000).
For glacial till sites, classical SPMs based on Richards' equation and the convection–dispersion equation (CDE) could somehow reproduce tile water outflow peaks but failed to describe the characteristic solute concentration pattern in the tile outflow associated with preferential flow (Wichtmann et al., 1998; Gerke and Köhne, 2004). In contrast, a 1D-DPERM representation of the field as a (vertical) soil block with lumped properties and assuming nonequilibrium-type preferential flow along soil macropores approximated both water and solute effluent patterns; here the observed bromide concentration peaks during leaching events could be described only when reducing mass transfer rate coefficient values (Villholth and Jensen, 1998; Gerke and Köhne, 2004). Based on the SPM concept, rapid appearance of Br in tile effluent was alternatively explained by lateral flow of perched water on a plow pan; this water eventually entered the subsoil layers in the disturbed soil regions close to tile lines (Kamra et al., 1999). Field tracer experiments analyzed with 2D-SPM and 2D-MIM showed that preferential flow as reflected by effluent concentrations from tile-drained field soils cannot be fully understood without considering 2D flow dynamics (Köhne and Gerke, 2005; Haws et al., 2005). A 2D model geometry representing a vertical cross-section perpendicular to the tiles may be the minimal requirement to describe soil water flow at a tile-drained field, which includes both vertical flow in the unsaturated zone and lateral flow toward tiles in the uppermost ground water zone. Moreover, a 2D-DPERM may allow considering spatial heterogeneity in both flow domains (Vogel et al., 2000). Still, the 2D-MIM could describe tile water flow and solute effluent concentrations in "weakly structured" soil at a marshland site (Köhne et al., 2006); however, the 2D MIM was not successful in simulating Br effluent concentrations at a glacial till site (Haws et al., 2005). By assuming two interacting flow domains, the 2D-DPERM could approximate water as well as pesticide leaching from tile outflow at a glacial till site (Gärdenäs et al., 2006). However, interpretation of preferential transport of adsorbing solutes is difficult without using conservative tracers to characterize flow (Flury, 1996). Thus, conservative tracers are also applied in combination with pesticides (Fox et al., 2006).
In most of the above-mentioned studies, model results were compared with the response at the tile-drain outlet point but not with measurements within the field. The simultaneous description of the flow process at the plot scale (i.e., tracer distribution in the soil profile) and the field-scale response (i.e., breakthrough curve in tile-drain effluent) requires distributed modeling or at least a procedure to distinguish between water and solute contributions from plot and field. If this is not warranted, the physical interpretation of such model results necessarily remains uncertain, particularly with model parameters obtained by calibration. Köhne and Gerke (2005) described a tracer experiment with Br strip application at the Bokhorst site, which included additional soil data and final tracer distributions in soil and plot-scale irrigation within the tile-drained catchment area. Here, a 2D-MIM gave a clearly better simulation of the observed preferential Br leaching than a 2D-SPM. However, explanations for soil water and tracer dynamics remained to some extent uncertain since a rather simplistic procedure was used for assessing plot- and field-scale contributions to tile outflow, and since inverse techniques were used for estimating water and transfer rate coefficients. Parameter fitting may lead to biased judgment on model performance since fitted "effective" parameter values may to some extent compensate for invalid model assumptions. Furthermore, hydraulic parameter values comparable to those of 1D-DPERM could not be used, and surface boundary conditions were simplified in the MIM simulations (Köhne and Gerke, 2005). In particular, the definition of the surface boundary conditions in the DPM allows distinguishing fractions of water and solutes that enter either the soil matrix domain or the preferential flow domain. Moreover, effects of using more realistic rainfall intensities instead of daily averages on preferential flow have not fully been explored.
The objective of this paper is to improve understanding of preferential flow in structured soils by analyzing the tile outflow and leaching patterns of the Br field tracer experiment (Köhne and Gerke, 2005) using a 2D-DPERM and a standard 2D-SPM for comparison without accounting for parameter fitting. We wanted to explore the dual-domain interactions in 2D and integrate 2D effects of tile drain, boundary conditions, and soil layering type of heterogeneity in modeling tracer leaching. The analysis distinguishes between effects of (i) model choice (2D-SPM vs. 2D-DPERM), (ii) surface boundary conditions for solute, (iii) solute mass transfer, and (iv) the contribution of plot-scale irrigation to field-integrated tile outflow and leaching using realistic temporal variation in irrigation intensities.
 |
Materials and Methods
|
|---|
Field Site and 1996–1997 Bromide Tracer Experiment
The Bokhorst tile-drained agricultural field site is located about 25 km south of the city of Kiel in northern Germany. Within this field, an area of about 5000 m2 is drained by a system of plastic (yellow PVC) drainpipes installed at about 1 m depth. The drains outlet flows were collected in a monitoring station equipped with a Venturi flume and an automated water sampler for determining outflow rates and bromide concentrations of the effluent. The Bokhorst site is gently sloped and had been used for a number of field tracer experiments (e.g., Wichtmann et al., 1998). For the 1996–1997 experiments analyzed here, a plot 10 m long (parallel to the drains) and 11.8 m wide (between two drains) was selected for tracer application and irrigation (Fig. 1
). At that plot, the surface was mostly flat with a slope of about 1%; the plot was located midway on a gently rising (1–2%) smooth hill slope. Instrumentation included a weather station (precipitation, air temperature and humidity recorded at 10-min intervals), piezometers, and tensiometers installed at 1, 3, and 5 m distance from the tile line in 0.15, 0.35, 0.5, 0.6, and 0.8 m depths. Winter wheat (Triticum aestivum L.) was growing on the field during the experiment. The heterogeneous soils originated from glacial till sediments, with a poorly drained Dystric Gleysol (FAO, 1998) as the most abundant soil type. The cultivated Ap horizon (0–0.3 m soil depth) had about 20% clay and 30% silt. Clay contents generally decreased with soil depth to about 10% at the 0.7- to 1.1-m depth. Soil structure ranged from subangular (0–0.3 m), angular blocky (0.3–0.95 m), and coherent (below 0.95 m depth). A plow-compaction pan was occasionally observed at 0.3 to 0.35 m depth. The hydraulic properties of the pan had been studied elsewhere (Köhne, 1999). However, its effects were not considered in the modeling here since the pan was only fragmentary and lacked continuity, as was observed on soil trenches north and south from the application strip. More details regarding the field site can be found in Lennartz et al. (1999) and Gerke and Köhne (2004). In the 1996–1997 experiment (Köhne and Gerke, 2005), the tracer was applied on a 10-m-long and 0.3-m-wide strip. The application strip was located at a horizontal distance of 0.85 to 1.15 m parallel to a tile line (Fig. 1 and 2
). A first Br application was on 18 Dec. 1996 (Day 0) when 3 L (i.e., 1 mm or 1 L m–2) of a solution containing 333 g L–1 Br (1.5 kg KBr dissolved in 3 L of water) was sprayed on the strip. Immediately after Br application, the plot area was irrigated with a backpack sprayer at 2 mm h–1 for 4 h. At that time, the wheat was in an initial phase (growth stage 10 according to Eucarpia (The European Association for Research in Plant Breeding) based on the scale of Zadoks et al. (1974). A second Br application followed on 25 Mar. 1997 (Day 97). This time, 6.4 kg of KBr salt (4.3 kg Br) was applied in solid form to the soil surface of the same 0.3-m-wide and 10-m-long strip. The plot area between two tile lines (Fig. 1) was then irrigated at a rate of 7 mm h–1 for 4.5 h in alternating 30-min irrigation intervals as follows. The one subplot of 10 by 5.9 m that included the Br application strip (bromide affected plot in Fig. 1) was irrigated five times (for 2.5 h total) and the adjacent 5.9 by 10 m subplot four times (for 2 h total). After a 12-h break, alternating irrigation at 7 mm h–1 continued in 30-min intervals for another 3.5 h on the following day (Day 98). The KBr salt had completely dissolved by the end of the 4.5-h irrigation on Day 97. One month earlier, on 24 Feb. 1997 (Day 68), the "initial" 2D distribution of Br concentration in the soil solution was determined from a 1.1-m-wide and 1-m-deep vertical cross-sectional area perpendicular to the tile line at the lower end of the application strip. At the end of the experiment, on 7 May 1997 (Day 140), the "final" 2D distribution of resident Br concentrations in the solution extracted from soil samples was again determined from a trench wall, this time at the upper end of the application strip (see Köhne and Gerke, 2005, for further experimental details).

View larger version (27K):
[in this window]
[in a new window]
|
FIG. 1. . The Bokhorst experimental site as used in the 1996–1997 tracer experiment, including the irrigated plot of 120 m2 and the application strip. The drainpipes are connected to the monitoring station (thick lines with arrows indicating the flow direction). The assumed boundary of the 5000-m2 catchment area of the drain outflow is shown as a hatched line.
|
|

View larger version (67K):
[in this window]
[in a new window]
|
FIG. 2. . Schematic picture of the two-dimensional (2D) vertical cross-section used in the simulations. Top: 2D single-porosity model simulation results for the pressure head distribution at Day 98.0 (i.e., between the two irrigation events). The colors represent soil matrix potentials, with pink indicating full saturation, and the arrows symbolize flow velocity vectors (colors representing Darcian flux in cm d–1) that are mostly small except in the vicinity of the drainpipe (the latter is indicated by a circle). Bottom: Finite element numerical grid.
|
|
The intention for strip Br application and plot irrigation within the drain catchment was to create conditions for simplifying the description of flow and tracer movement in a 2D vertical cross-section perpendicular to the tile lines. In contrast to the aerial Br application of previous study (Wichtmann et al., 1998), the strip experiment allowed using 2D soil information and residual Br spatial distribution along soil trenches. Furthermore, we intended to assess and better explain the data by using both the information at the plot scale and the response at the tile outlet representing the field scale. Splitting the plot in two areas of alternating irrigation was done for technical reasons; we could only use an irrigation device of maximal 6 m width.
Two-Dimensional Numerical Analysis: Single- and Dual-Permeability Models
The SWMS II finite element numerical code (Vogel, 1987) was used as the reference conventional 2D-SPM. SWMS II solves the Richards' equation for describing 2D Darcian water flow and the advection–dispersion equation for calculating solute movement.
The 2D-DPERM numerical model applied here was first described in Vogel et al. (2000) and is based on the concept of Gerke and van Genuchten (1993a) for water and solute movement in a variably saturated structured soil. The DPERM assumes that a structured porous medium can conceptually be separated into two mobile pore domains having two sets of water retention and hydraulic conductivity functions, one for the matrix (subscript m) and another for the fracture (subscript f) pore system (here denoted as preferential flow [PF] domain). The "local" (i.e., pore domain) properties are related to those of the total soil volume by
 | [1a] |
 | [1b] |
 | [1c] |
 | [1d] |
 | [1e] |
where wf is the relative volumetric proportion of the preferential flow domain, wm = 1 – wf,
is porosity (L3 L–3),
is volumetric water content (L3 L–3), q is fluid flux density vector (L T–1), c is solute concentration (M L–3), and Ks is saturated hydraulic conductivity (L T–1). Note that Eq. (1e) assumes equilibrium saturation in both domains (i.e., hs,m = hs,f = 0). The 2D flow of water in the dual-porosity medium is described with two coupled Richards' equations:
 | [2] |
where h is the pressure head (L); K is hydraulic conductivity tensor (L T–1); C is the specific water capacity (L–1); S is a sink term for root water uptake (T–1), here assuming that roots extract water only from the matrix domain; z the vertical coordinate, taken to be positive upward (L); t is time (T); and
w is the transfer term (T–1) for water exchange between the fracture and the matrix pore systems assumed (Gerke and van Genuchten, 1993a) as
 | [3] |
where
w is the first-order water transfer rate coefficient (L–1 T–1), here taken (Ray et al., 2004) as
 | [4] |
where
ws is the water transfer rate coefficient at saturation (L–1 T–1), and Kar is a relative unsaturated hydraulic conductivity function (–) of the interface (i.e., the aggregate surface, subscript a) between the soil matrix (SM) and the PF domain. Values of Kar range from 0 to 1 taken as the minimum of the PF and SM domain conductivities evaluated for upstream pressure (i.e., Kar = min[Kfr(hf)],[Kmr(hf)] for hf
hm and Kar = min[Kfr(hm)],[Kmr(hm)] for hf < hm). Values of
ws (Eq. [4]) can account for mass transfer reductions due to, for instance, aggregate coating effects (Gerke and Köhne, 2002) in a more lumped manner as in a previously suggested formulation (Gerke and van Genuchten, 1993b):
 | [5] |
where a is the characteristic radius or half-width of the matrix structure (L), ß is the dimensionless geometry coefficient, Ka is hydraulic conductivity depending on the degree of SM-domain water saturation (L T–1), and
w = 0.4.
The sink term, S (Eq. [2]), describes root water uptake according to Feddes et al. (1978). We assumed a vertically linear root distribution with maximum at the soil surface and minimum at a depth of 1 m. Relatively low sensitive parameter values of the plant water stress function (Feddes et al., 1978) were used (h1 = +10 cm, h2 = +9 cm, h3 [high] = –2000 cm, h3 [low] = –8000 cm, and h4 = –120,000 cm). Evapotranspiration (ET) rates were estimated according to the water budget using potential ET rates (Haude, 1955) obtained from field data of temperature and relative humidity (Köhne and Gerke, 2005, p. 82). Relatively low actual evapotranspiration rates (3–45 mm mo–1) were obtained. There was a period with frost in December and January, however, with some snow forming a discontinuous and negligibly thin snow cover. Thus, no plant water stress occurred during most of the experimental period covering cool and humid winter and early spring months.
The 2D transport of solutes including linear equilibrium sorption in a dual-permeability soil involves two coupled advection–dispersion equations (e.g., Vogel et al., 2000) as
 | [6] |
where D is the dispersion coefficient tensor (L2 T–1), R is the dimensionless retardation factor (here, R = 1), Sc is a sink term for root uptake of solutes (M L–3 T–1) describing passive Br uptake with water (i.e., Sm·cm). Note that root uptake of solutes was not further considered here because it had only negligibly small effects on final Br distribution, effluent concentrations, and mass balance. The term
s in Eq. [6] is the solute mass transfer term (M L–3 T–1) evaluated as (Gerke and van Genuchten, 1996)
 | [7] |
in which ci = cf if water is transferred from the PF into the SM domain (i.e.,
w positive) and ci = cm if water transfer is in the opposite direction,
s =
s'wm
m or
s = wm
m(ß/a2)Da is the first-order solute mass transfer rate coefficient (T–1), here expressed comparable to Eq. [4] as
 | [8] |
where
ss is the solute transfer rate coefficient at saturation (T–1). In this work, the relative saturation
ar =
a/
as of the interface between the SM and the PF domain is assumed to be equal to the (more dynamic) relative saturation of the PF domain (Ray et al., 2004). In analogy to Eq. [5],
s' is of the form (van Genuchten and Dalton, 1986):
 | [9] |
where ß and a are the same as for water (Eq. [5]), and Da is the effective diffusion coefficient (L2 T–1) of the exchange term.
The soil hydraulic functions were described using the modified version of the van Genuchten–Mualem formulation (van Genuchten, 1980) to avoid numerical problems of the water contents near saturation as (Vogel and Cislerova, 1988; Vogel et al., 2001)
 | [10] |
 | [11] |
where
 | [12] |
with
 | [13] |
where
s are saturated and
r residual water content parameters (L3 L–3), respectively; n (–) and
VG (L–1) are empirical shape parameters with n > 1 and m = 1 – 1/n; Se = (
–
r)/(
s –
r) is effective saturation (–), and Ks is saturated (L T–1) and Kr relative hydraulic conductivity (–). In this modified expression,
s is replaced by a fictitious parameter,
s*, with the value of
s* being slightly higher than that of
s. Parameter hs denotes the minimum capillary height (L) at
s*. In this study, we used the same soil hydraulic parameters (Table 1) as in previous 1D dual-permeability simulations (Gerke and Köhne, 2004), which were derived from calibration of an earlier tracer experiment (Wichtmann et al., 1998). The composite parameters for the 2D-SPM were obtained using Eq. [1b–1e]. For the subsoil (40–200 cm depth), an anisotropy of hydraulic conductivity of Kxx/Kzz = 3 was assumed according to lab-determined values for vertical and horizontal hydraulic conductivities (Köhne, 1999, p. 133).
View this table:
[in this window]
[in a new window]
|
TABLE 1. Hydraulic parameters used in the simulations of the Bokhorst bromide field tracer experiment for the single-domain model (2D-SPM) and the dual-permeability (2D-DPERM) model for soil matrix (SM) and preferential flow (PF) domains assuming wf = 0.05 for both layers. Anisotropy of hydraulic conductivity (Kxx/Kzz = 3) was assumed in 40–200 cm depth. Solute transport parameters were fixed as L = 20 cm, T = 2 cm, and D0 = 1.2 cm2 d–1.
|
|
The solute transport parameter (dispersion coefficient tensor D [L2 T–1]) was evaluated (Bear 1972) as
 | [14] |
where Do is the molecular diffusion coefficient (L2 T–1),
L is the longitudinal and
T the transversal dispersivity (L), q is the magnitude of the vector of water flux (L T–1),
ij is the Kronecker delta function (
ij = 1 if i = j, and
ij = 0 if i
j), and
=
7/3/
s2 (Millington and Quirk, 1961) is a dimensionless tortuosity factor. Here, transport parameters
L = 20 cm,
T = 2 cm, and D0 = 1.2 cm2 d–1 were used for both layers and domains without further analyses.
The transfer term parameters (Table 2) determining the exchange of water and solute between the domains (Eq. [4] and [8]) were assumed here to be
ws = 0.01 and
ss = 0.2 (0.02) for the Ap horizon and
ws = 0.005 and
ss = 0.1 (0.01) for the CSd horizon and the subsoil. These values are slightly higher as those of the 1D-DPERM analysis (Gerke and Köhne, 2004), where values of
ws =
ss = 0.006 for the Ap and of
ws =
ss = 0.004 for the CSd horizon were used. The lumped coefficients of the 1D analysis were based on values of
w = 0.4 (Gerke and van Genuchten, 1993b), ß = 15 (Ap) and ß = 10 (CSd), and a = 1 cm in Eq. [5] and [9] and assuming 1000 times reduced hydraulic and diffusive exchange coefficients compared with the SM domain parameters. The values for the coefficients representing the soil structural geometry, ß and a, were selected according to qualitative descriptions of the soil structure (Gerke and Köhne, 2004). The smaller reductions in hydraulic and diffusive exchange rate coefficients (less than 100 times) compared with 1D-DPERM analyses (1000 times) were suggested from results of studies with intact and removed coatings (Gerke and Köhne, 2002; Köhne et al., 2002a). The larger values for solute exchange rate coefficient were similar to those of the MIM simulations (Köhne and Gerke, 2005); a second set of smaller
ss values was tested since exchange of solutes was found to be more sensitive than that of water for the Bokhorst site.
View this table:
[in this window]
[in a new window]
|
TABLE 2. Dual-permeability (DPERM) mass transfer rate coefficients (relative) for water exchange, ws, and for solute exchange, ss, used in the simulation scenarios for both alternative surface boundary conditions: Br addition to the soil matrix domain only and addition to both pore domains.
|
|
Equations [2] and [6] were solved numerically using fully implicit Galerkin finite elements. Discretization of the 2D vertical cross-sectional domain of x = 5.9 m width (i.e., half-distance between tiles) by z = 2 m depth was by 4464 nodes and 4331 rectangular elements (Fig. 2). Smaller nodal distances (0.0025–0.02 m) were selected near the soil surface and around the tile drain. Nodal increments elsewhere in the model domain increased up to 0.2 m. Mass balance errors were always smaller than 1% for both water and solutes. Flow and transport equations were solved sequentially for the two domains, beginning with the PF domain and then solving for the SM domain. The transfer terms were evaluated at the "j-1 time level" where j is the "present time level." As for 1D simulations, higher water transfer rates led to model failure.
Initial and Boundary Conditions and Scenarios
For water flow, we used a pressure head profile according to measured tensiometer values as initial condition on Day 0. For the runs starting on Day 68, the initial pressure head distribution was derived from simulations performed for the period Day 0 to Day 68. The initial condition for Br transport was a solute-free 2D cross-sectional soil profile for Day 0. The spatially interpolated 2D distribution of measured Br concentrations was used for simulations starting at Day 68.
At the soil surface, atmospheric boundary condition (BC), flux type, was imposed to represent rainfall- and irrigation-induced infiltration and evapotranspiration (Köhne and Gerke, 2005). For conditions when rainfall or irrigation rates exceeded the maximal infiltrability of soil at the surface, which happened twice during the irrigation periods, on Days 97 and 98, surplus water was temporarily stored and allowed to resume infiltrating after the end of the irrigation event assuming a flux-type condition. For 2D-DPERM, the real 0.5-h switching intensities were used to simulate irrigation. By contrast, for 2D-SPM, intensities averaged over 4.5 h and 3.5 h irrigation durations, respectively, were used since the higher real intensities would have resulted in significant ponding, creating water mass balancing problems. The effect of using averaged irrigation intensities was negligible for 2D-SPM solute transport simulations.
At the Br application strip, a third-type solute BC with prescribed Br concentrations of the infiltrating water was imposed. The Br concentration of the application was calculated using the cumulative infiltration flux and total applied Br mass only for the irrigation during Day 97. The total load of Br was identical for all simulation cases and both BC scenarios (i.e., Br applied [i] to both the PF and SM domains and [ii] to the SM domain only). In 2D-DPERM simulations, the lower SM hydraulic conductivity resulted in relatively small cumulative irrigation fluxes to the SM domain and water redistribution at the surface from the SM to the PF domain. In that case, the Br concentration of the irrigation water entering the SM domain was increased to match the balance load (4.3 kg per 10-m-long application strip). For the "SM domain only" scenario, Br was exclusively entering the SM domain with the irrigation water. For the "both domains" scenario, about 5 mm of irrigation water (out of 40 mm) were forced to infiltrate after the end of the irrigation event on Day 97, leading to an increase in influx concentration from a value of 42.8 mg cm–3 (initially calculated) to 51.3 mg cm–3 (required for balancing).
At the bottom boundary, a no-flow condition as in Köhne and Gerke (2005) was imposed in 2 m depth (Fig. 2). At the left and right sides of the domains, a no-flow condition was imposed assuming symmetry perpendicular to drain tiles. Tile-drain effluent of water was simulated as a seepage face boundary condition (i.e., Dirichlet pressure-type condition) with drainage beginning at saturation (h
0 cm) and Neumann flux-type condition with zero flux for unsaturated conditions that was imposed across the inner tile diameter of 0.04 m (Köhne and Gerke, 2005). Any possible drainage resistance by drainpipes was not explicitly considered. The tile was represented by means of three numerical nodes at 1 m depth, forming the effective length of the drain as 4 cm. For simulating Br– leaching out of the drain tile, a zero-gradient solute BC was ascribed for both the PF and the SM domains (i.e., Br– may enter the tile with the effluent water).
Scenario simulations of the total experimental period (0–140 d) were used to obtain initial conditions (i.e., starting with conditions on Day 68) for the specific numerical experiments that focused on the period between 2 Feb. 1997 (Day 69) and 7 May 1997 (Day 140). For the DPERM, relatively high and low solute mass transfer rate coefficients and two alternative upper boundary conditions for Br– application in both domains and only in the matrix domain (Table 2) were simulated. Simulations with both 2D-SPM and 2D-DPERM were restarted at Day 68 with measured Br– distribution as initial condition for solutes for comparing results of 2D-SPM and 2D-DPERM for the four cases (i.e., two boundary conditions, Br addition to SM domain and to both domains, for the respective two mass transfer rate coefficients; Table 2).
Plot Scale versus Field Scale
For comparing model results from 2D cross-sections with observed field-scale data, simulated tile outflow of water and leaching of Br as calculated for a 2D vertical cross-sectional area (Fig. 2) were assumed to be representative of one-half of the irrigated horizontal plot area (Fig. 1). The 2D simulations were performed separately for irrigated plots (i.e., by applying rainfall and irrigation rates) and for nonirrigated plots (i.e., by applying rainfall only), the latter to represent drainage and Br leaching from the rest of the tile-drained area. The two drainage time series from irrigated and nonirrigated areas were then combined to obtain the area-weighed tile outflow and Br effluent concentration for comparison with the measured data. The contribution of the irrigated plot of 120 (60) m2 to drainage (leaching) from the total field of about 5000 m2 (Fig. 1) was considered by means of several weighing procedures, including weighing from dual-porosity concept, from the 2D concept and the partial-area irrigation and strip application of Br (see Appendix). All simulation results for the 2D-SPM and 2D-DPERM are presented as field-scale area weighed. The "Nash criterion" describing the model efficiency coefficient, E, was calculated (Nash and Sutcliffe, 1970) as
 | [15] |
where Xo is observed variable, Xm is modeled variable and is the mean value of observed variable.
 |
Results
|
|---|
For the simulated 2D vertical cross-section of 2 m depth and 5.9 m length (Fig. 2), the 2D-SPM results for flow on Day 98.0 (i.e., between the two irrigation periods) reflect a hydraulic situation for the near-saturated topsoil and the fully saturated subsoil with flows directed toward the tile, although drainage rates are still relatively small. Flow in the topsoil is directed mostly laterally here (Fig. 2, top), instead of an expected vertical flow direction in unsaturated soil. This is because relatively small flow intensities (horizontal
0.4 cm d–1 and vertical
0.01 cm d–1) were computed for the shallow fringe near saturation just at the boundary between the Ap and CSd horizons on top of a gently sloped water table; moreover, no rainfall occurred at this time. With increasing thickness of the shallow unsaturated zone in the vicinity of the drain, the flow direction becomes more vertically oriented. Lateral flow dominates in the saturated zone, where flow rates increase toward the drainpipe and in the region the water table below the Br application strip. Drain water outflow (Fig. 3
) could be described similarly well with both the 2D-SPM and the 2D-DPERM, using the parameter setting of Gerke and Köhne (2004). The values of the Nash–Sutcliffe criterion (Table 3) for drain outflow were E = 0.4 (2D-SPM) and E = 0.598 (2D-DPERM) for the total period (0–140 d).

View larger version (27K):
[in this window]
[in a new window]
|
FIG. 3. . Time series of tile outflow comparing data and model results for a two-dimensional single-porosity model (2D-SPM) and two-dimensional dual-permeability model (2D–DPERM). Drain discharge rates (lines) as well as daily rainfall and plot scale irrigation rates (blue columns at the top) are shown (a) for the total experimental period, and in greater detail (b) for the first and (c) for the second drainage period including the irrigation (Days 97 and 98) and a larger rainfall event (Day 100); rates for (c) are on a 12-h basis.
|
|
View this table:
[in this window]
[in a new window]
|
TABLE 3. Values of the Nash–Sutcliffe model efficiency coefficients, E, comparing drain discharge flow and cumulative Br leaching curves with results of the simulation scenarios.
|
|
Both approaches could reproduce especially the larger drainage peaks indicated by increased values of the Nash criterion of E = 0.69 (2D-SPM) and E = 0.71 (2D-DPERM) for the rainy period (Days 55–80). The models overestimated drainage rates during some smaller precipitation events during winter (Fig. 3a, Days 20–50) and later in spring (Fig. 3a, Days 90, 110). The more detailed picture for rainfall period between Days 50 and 85 (Fig. 3b) shows that both 2D-SPM and 2D-DPERM give an excellent description of the data. Both models approximated drain discharge during the irrigation period between Day 97 and Day 98 (Fig. 3c) when using area-weighed field-scale values (see Appendix). Correspondence is higher for 2D-DPERM (i.e., E = –0.05 [2D-SPM] and E = 0.55 [2D-DPERM] for Days 97–103; Table 3). The match between simulated and measured data could most likely be improved by manipulating model parameters. However, model calibration was beyond the scope of this study. Note that the daily irrigation rates (circled in Fig. 3c) are presented as the plot scale values. In contrast to the relatively small discharge rates that were observed following local irrigation of the plot of 120 m2, the considerably larger drain discharge after Day 100 (Fig. 3c) resulted from field-scale rainfall (i.e., related to the total drain catchment area of 5000 m2). The plot scale effect on tracer leaching is discussed below.
Results of the 2D-DPERM simulations of pressure heads, water fluxes, and water transfer rate coefficients (Fig. 4
) are shown as vertical cross-sectional images directly after the first irrigation (at Day 97.74). At this time, the PF domain was almost entirely saturated compared with the SM domain, which was only saturated near the surface and in the subsoil (Fig. 4). Flow velocity vectors represent local water flux within the PF domain. Here, the scale of Darcian flow velocities (Fig. 4, top right) is about 40 times larger than the one used for 2D-SPM simulations (Fig. 2). However, because of the small volume fraction of the PF domain, the simulated specific discharge is comparable in both SPM and DPERM (Fig. 3). Differences between the pressure head distributions of the SM and PF domains (Fig. 4, top) indicate spatially varying local nonequilibrium conditions. All along the cross-section in the upper part of the subsoil (Fig. 4, bottom), water is transferred from the more saturated PF domain to the SM domain (i.e., positive mass transfer rates). The green-colored regions (Fig. 4, bottom) indicate conditions near pressure head equilibrium with negligibly small water transfer. Positive exchange rates (blue color) indicate transfer from the PF to the SM domain. The horizontal variation in the 2D distribution of water transfer rates here solely results from effects of the domain-specific flow fields close to the tile since spatial variability of hydraulic properties was not further considered.

View larger version (23K):
[in this window]
[in a new window]
|
FIG. 4. . Vertical cross-sectional images of water flow simulation results obtained with a two-dimensional dual-permeability model (2D-DPERM) for a nonequilibrium situation at Day 97.74 (i.e., directly after the first irrigation) showing (top) pressure head distribution in the soil matrix (SM) (top left) and preferential flow (PF) domain (top right) with velocity vectors, and (bottom) 2D water transfer rate distribution. Darcian flow velocities in cm d–1 are indicated by the color of the arrows.
|
|
The tensiometer data (Fig. 5
) observed at 35 and 50 cm soil depth and 1 m distance from the tiles tended toward positive pressure heads during first irrigation (Day 97) and the rainfall event (Day 100). Data during the second irrigation were not available due to technical failure. Tensiometer values from other depths and distances from the tile line (Köhne and Gerke, 2005) show similar curves. The tensiometer data from two depths indicate that the water table was below the plow pan in about 40 cm soil depth during first irrigation; during the following rain event, a second perched water table may have occurred above the plow pan (if present at all at about 30 cm depth). The perched water may have either stimulated water entering larger pores of the PF domain in the pan or initiated lateral water movement on the pan. The real influence of perched water on a plow pan cannot be interpreted from these tensiometer data; only the 2D vertical distribution of residual Br concentrations may indicate some effect (see below). In the field, neither a perched water table nor any ponding on a plow pan could be observed (piezometers were not installed in such shallow depths).

View larger version (30K):
[in this window]
[in a new window]
|
FIG. 5. . Time course of pressure heads at 35 cm (top) and 50 cm (bottom) depth and 1 m apart from the tile line, comparing measured data with simulated values obtained with two-dimensional single-porosity model (2D-SPM) and the two-dimensional dual-permeability model (2D-DPERM) for soil matrix (SM) and preferential flow (PF) domains. The vertical lines at Day 97.74 indicate the time directly after first irrigation, which is comparable to the 2D cross-sections.
|
|
In contrast to the data, the pressure head values simulated with 2D-SPM increased and decreased more sharply and reached higher values during irrigation but lower values during the rainfall event (Fig. 5). The pressure heads obtained with 2D-DPERM for the PF domain are somewhat similar to those of 2D-SPM; the values for the SM domain are more similar in shape to those of the data, although at a significantly lower level. For both soil depths, pressure head values increase faster and higher in the PF domain compared with the SM domain during infiltration events. Despite its relatively high permeability, the restricted PF domain pore fraction receives water from upper soil layers at rates still higher than downward flow and transfer into the SM domain. Saturation differences between the PF and SM domains and nonequilibrium conditions are larger in the greater depth during irrigation (Fig. 5, bottom). The 2D-DPERM predicted small water saturations for the second irrigation event and only for the PF domain; the SM domain remained largely unsaturated. Note that no parameter calibration was performed to improve fit between measured and simulated values. The 2D-DPERM predicts local nonequilibrium in pressure heads between the SM and PF domains for the total period (Fig. 5), while the direction of mass transfer is changing during and after irrigations. Before the irrigation, the less-negative pressure heads in the SM domain (Day 97, Fig. 5) result from the relatively faster draining of the PF domain. Consequently, before irrigation, water transfer at the soil surface is directed from the SM domain into the less-saturated PF domain, while during irrigation (or always at greater depths), transfer is from the PF to the SM domain. Furthermore, water transfer from the PF into the SM domain is restricting the increase in the PF domain's pressure head such that maximal values remain smaller, as those of the 2D-SPM simulations (Fig. 5).
The 2D water saturation distribution at Day 97.74 for the cross-section (Fig. 4) corresponds to the one at the two depth indicated by lines in Fig. 5. At this time, simulated domain-specific pressure heads are close to local equilibrium in 35 cm depth and nearly identical with data, and still unsaturated. For the 50 cm depth, the 2D-DPERM predicts local nonequilibrium while the tensiometer data probably tend to reflect mainly the pressure in the larger pores.
Results of the 2D-DPERM solute transport simulation at Day 97.74 shortly after the first irrigation show a situation of local nonequilibrium in Br concentrations (Fig. 6
). At the soil surface, all Br was applied to the SM domain, in which transport and dispersive spreading was relatively small (Fig. 6, top left). In the PF domain, the Br plume extended to greater depths, where the shift in the plume indicates the beginning of lateral Br movement toward the tile near the water table (Fig. 6, top right). The Br displacement in the PF domain generally reflects the corresponding preferential water flow pattern (Fig. 4), modified by dispersion and solute mass transfer. In the upper part of the topsoil, solute mass transfer (Fig. 6, bottom) directed from the SM into the PF domain is relatively high, as indicated by negative rates (green and red colors). Here, the "diffusive" mass transfer term component (Eq. [7]) is dominating, while the "convective" water transfer–related component is negligible (compare Fig. 4). Once Br entered the PF domain, it was rapidly transported vertically downward such that enhanced Br leaching in tile effluent originating from application on Day 97 just started around this time (Day 97.75; see below). In the subsoil, solute transfer is directed from the PF into the SM domain (Fig. 6, bottom), as indicated by positive rates (pink color). The blue-colored regions indicate conditions near equilibrium in concentrations with negligibly small solute transfer. Although transfer proceeds at smaller rates than in the topsoil, it nevertheless affects the shape of the Br plume in the PF domain. The sharp front in positive solute mass transfer rates corresponds with the layer boundary at 30 cm soil depth and is indirectly caused by the corresponding discontinuity in water contents.

View larger version (36K):
[in this window]
[in a new window]
|
FIG. 6. . Vertical cross-sectional images (2-m-wide section below application strip) presenting results of the two-dimensional dual-permeability model (2D-DPERM) solute transport simulation during nonequilibrium situation at Day 97.74 (i.e., shortly after the first irrigation) for (top) 2D distribution of Br concentrations in (top left) soil matrix (SM) and (top right) preferential flow (PF) domain, and (bottom) for 2D vertical distribution of solute mass transfer rates. Here Br was applied in the SM domain only, and solute transfer rate coefficient, ss, was assumed relatively high (Table 2).
|
|
We compared measured Br leaching in tile-drain effluent with simulated results obtained with 2D-SPM and 2D-DPERM for Br application in both domains (Fig. 7
) and Br application in the SM domain only (Fig. 8
). Here, simulations started at Day 68, assuming a linearly interpolated distribution of measured Br concentration (c.f. Fig. 6 in Köhne and Gerke, 2005); nodes below 1 m depth received a zero initial Br concentration condition. The 2D-SPM overestimates the Br effluent concentration (Fig. 7 and 8) after the first irrigation. Later, after the second irrigation, effluent concentrations are underestimated such that cumulative leached mass of Br as calculated with 2D-SPM is significantly lower than as measured. Although discharge rates were overestimated, practically Br-free water from other regions of the cross-section contributed to the diluted Br concentration pattern in case of 2D-SPM. Furthermore, Br concentrations were simulated only for a single 2D cross-section of one plot regarded as representative for the three-dimensional (3D) flow field that is affecting the tile effluent.

View larger version (25K):
[in this window]
[in a new window]
|
FIG. 7. . Time series comparing measured and simulated results of (top) Br concentrations in tile outflow and (bottom) cumulative Br mass leaching as in Fig. 7, but here for Br application in the soil matrix (SM) domain only. 2D-SPM, two-dimensional single porosity model; 2D-DPERM, two-dimensional dual-permeability model.
|
|

View larger version (26K):
[in this window]
[in a new window]
|
FIG. 8. . Time series comparing measured and simulated results of (top) Br concentrations in tile outflow and (bottom) cumulative Br mass leaching. Curves are obtained with two-dimensional single porosity model (2D-SPM) and with two-dimensional dual-permeability model (2D-DPERM) for higher (2D-DPERM high) and lower (2D-DPERM low) values of the solute mass transfer rate coefficient, ss, for the case of bromide application in both domains (Table 2). Blue columns indicate irrigation and a following rainfall event.
|
|
In contrast, the 2D-DPERM for the case of Br application in both domains (Fig. 7, top) largely overestimates Br concentrations in tile effluent, particularly after the first irrigation. Reduced mass transfer rate coefficients (Fig. 7, bottom) resulted in the largest Br leaching, mainly because solute transfer into the SM domain of the subsoil is restricted. On the other hand, effluent concentrations are reduced when solute transfer is increased, which results from stronger Br "absorption" by the SM domain of the subsoil. Still, 2D-DPERM simulations strongly overestimate Br leaching during the first event because in this scenario, a significant fraction of Br could directly enter the PF domain at the soil surface (Fig. 7). Note that the cumulative Br mass was assumed to start from zero at Day 95.0, not accounting for the previous Br leaching episodes.
For Br surface application to the SM domain only (Fig. 8), the effluent concentrations are significantly better approximated by the 2D-DPERM approaches. The best fit for solutes according to Nash criterion with E = 0.759 was found for 2D-DPERM for the scenario assuming smaller exchange coefficient and Br addition and surface infiltration into the SM domain only. The values of E in Table 3 may all appear relatively low; however, we did not intend here to further improve model performance by parameter fitting. A better fit and non-negative values of E for Br would be obtained simply by assuming some Br background concentration.
For both the relatively high and low solute transfer rate coefficients, the characteristic feature of a higher Br concentration during the second irrigation is now captured. In contrast, the 2D-SPM results (i.e., the same as in Fig. 7, but different scale) fall short for the second irrigation and following periods. The 2D-DPERM results could possibly be further improved; however, no parameter calibration was attempted here. Note that the relatively small Br concentration peak in tile effluent at Day 100 (Fig. 7 and 8) follows the field-scale rainfall event (c.f. Fig. 3). Especially for the second irrigation event, the 2D-DPERM approach (Br application in matrix domain and high solute transfer) approximates Br concentrations significantly better than 2D-SPM (Fig. 8, top). Cumulative Br leaching was in the same order of magnitude as the measured values; the slight overestimation of the leached mass (Fig. 8, bottom) is caused partly by slightly overestimated water effluent rates (Fig. 3). The best match of Br concentrations according to Nash criterion with E = 0.66 was found for 2D-DPERM (low, matrix) assuming soil surface infiltration in the SM domain only. While some of the values of E (Table 3) may appear quite low, we did not intend to further calibrate model parameters or boundary conditions. Moreover, model calibration would have required to consider Br background concentrations, which were not exactly known.
As expected, the simulated 2D vertical distribution of residual Br concentrations in the profile (Fig. 9
) at the end of the experiment (Day 140) could not fully resemble the heterogeneous pattern of measured Br concentrations since the models did not include spatially distributed hydraulic parameters, among other factors. Nevertheless, we may compare the location of the peak Br concentration and the overall shape and spreading of the Br plume for model evaluation. Compared with the data (Fig. 9a), the 2D-SPM (Fig. 9b) predicts a deeper leaching of the Br plume and a similar spreading at generally lower values of the Br concentration in the center of the plume. In contrast, the 2D-DPERM for the case of Br application to the SM domain and low mass transfer predicts a Br plume peak (Fig. 9e) that is not as deep and shows less spreading but higher Br peak concentrations, indicating a somewhat better approximation toward the data compared with 2D-SPM. The Br concentration plume for the SM domain of the 2D-DPERM (Fig. 9c) mainly resembles that of the bulk soil volume (Fig. 9e); the distribution of Br concentrations in the PF domain (Fig. 9d) reflects rapid vertical transport and lateral-oriented movement toward the drain tile in the subsoil. Although somewhat similar in shape, the 2D-DPERM-plume (Fig. 9e) is significantly deeper and wider than the one obtained with the MIM approach (c.f. Fig. 9 in Köhne and Gerke, 2005). Furthermore, we note that the measured data (Fig. 9a) indicate a sharp contrast in Br concentrations across the interface between plow layer and subsoil, which apparently was not found in the 2D-SPM (Fig. 9b) and 2D-DPERM (Fig. 9c) simulations. However, more detailed interpretations of the differences between the simulated and observed Br plume, including possible effects of soil heterogeneity and of perched water table at a plow pan, would require further measurements and are beyond the scope of this study.

View larger version (33K):
[in this window]
[in a new window]
|
FIG. 9. . Two-dimensional spatial distribution of Br concentration in the soil solution at the end of the experiments on Day 140 (7 May 1997) of transect 2 for (a) observations (DATA), and simulations obtained with (b) the two-dimensional single porosity model (2D-SPM), (c) the two-dimensional dual-permeability model (2D-DPERM) for the soil matrix (SM) domain, (d) the 2D-DPERM for the preferential flow (PF) domain, and (e) the 2D-DPERM composite Br concentration (related to bulk soil volume, Eq. [1d]). The 2D-DPERM simulations assumed Br application in the SM domain only and relatively low mass transfer rate coefficient.
|
|
 |
Discussion
|
|---|
The results allow discussion of problems related to both the modeling aspects (basic assumptions, boundary conditions for solutes, 2D effect) and the specific experimental conditions (strip application, lateral plume, plot vs. field scale). Although we believe that water flow description is possible with both approaches, we do not intend to distinguish models on the drain discharge curves alone without considering tracer transport.
Effects of Model Approaches
In this study, we reanalyze a tracer experiment from a field by considering a model with two interacting flow domains (2D-DPERM) that allows for local nonequilibrium in both pressure heads and concentrations to consider soil structure effects. As a reference, a single porosity approach (2D-SPM) is used, which assumes local equilibrium conditions and "effective" parameters to include soil structure effects.
Additional comparisons with simulation results using previously studied 1D-DPERM and 2D-MIM approaches were not included here. For 1D-DPERM, the experimental design consisting of the strip Br application as a KBr salt (Köhne, 1999) did not justify a 1D vertical approach as for a previous Br tracer experiment (Wichtmann, 1994) where Br was applied across the whole catchment and used for the 1D-DPERM study (Gerke and Köhne, 2004). The 2D-MIM was not included in this study because the surface boundary conditions for particular scenarios were not compatible with those imposed in DPERM: Whereas in DPERM, solute can enter the SM domain (i.e., immobile region in MIM), all solute can only enter the mobile region (i.e., PF domain) in MIM.
However, at the plot scale, 1D and 2D cases for SPM and DPERM can be compared for the case of uniform Br application along the soil surface. The scenario (Fig. 10
) assumes Br addition to the SM domain, larger mass transfer rate coefficients, and a zero initial Br concentration. For water, the tile-drain outflow curves demonstrate largest peaks following infiltration for both 1D compared with 2D approaches, but also larger discharge peaks for both SPM approaches compared with the respective 1D- and 2D-DPERM approaches (Fig. 10, top). These characteristic differences result from additional storage effects in 2D due to elevated and curved water tables caused by tile drainage not considered in 1D for the former and from water exchange effects in the DPERM for the latter. The Br concentrations in the drain effluent predicted by both 1D- and 2D-DPERM models are initially larger compared with both SPM approaches (Fig. 10, bottom). Later, Br concentrations for 2D-DPERM remain below 1D-DPERM because of mixing and dilution during lateral movement in the case of 2D. In contrast, the Br effluent concentrations are initially larger for 2D-SPM compared with 1D-SPM. The results (Fig. 10) reflect effects of lateral transport toward the drain in 2D that lead to higher Br concentrations during drainage events in 2D-SPM when lateral flow increases within an elevated groundwater zone.

View larger version (23K):
[in this window]
[in a new window]
|
FIG. 10. . Time series comparing simulated plot scale results of (top) drain discharge and (bottom) Br resident concentration at the drain position obtained with single (one-dimensional single porosity model [1D-SPM] and two-dimensional single porosity model [2D-SPM]) and dual-permeability models (one-dimensional dual-permeability model [1D-DPERM] and two-dimensional dual-permeability model [2D-DPERM]), assuming zero Br concentration at Day 97 (higher mass transfer rate coefficient, ss, and Br addition to the soil matrix domain only for DPERM).
|
|
As a continuation of 1D-DPERM modeling, the same parameter values were used in the 2D-PERM approach here. However, these values cannot be used in the MIM simulations. Previous work suggested that MIM could be used after fitting exchange coefficients (Köhne and Gerke, 2005). However, parameter fitting may lead to biased judgment on model performance, since fitted "effective" parameter values may to some extent compensate for wrong model assumptions, analogous to the dispersion length as an "effective" parameter, which in 1D-SPM is increasing with distance and water saturation if physical nonequilibrium conditions exist in fine-textured soil (Vanderborght and Vereecken, 2007). Here, for better process understanding, a different concept is tested based on the physical consideration that for the particular soil structures, macroscopic water movement occurs in the SM domain. Matrix involvement in the overall flow process was suggested based on observations of soil wetting during tension infiltrometer experiments and Brilliant Blue dye tracer tests, which at this site showed transport in the soil matrix, particularly in the top soil horizon. Dye tracer tests usually result in dye transport constricted largely to preferential flow paths, due to adsorption and retardation of dye in the soil matrix. Therefore, dye tracer tests normally may lead to underestimation of "vertical" conservative solute transport in the SM domain (Kasteel et al., 2005). This underlines the need to consider matrix tracer transport at the Bokhorst site.
The 2D-MIM (Köhne and Gerke, 2005) did not yield closer correspondence with data as 2D-DPERM shown here, despite parameter optimization. In particular, the overestimation of the first concentration peak during the irrigation period by 2D-SPM (Fig. 8) and 2D-MIM (Fig. 7d in Köhne and Gerke, 2005) may similarly indicate a conceptual limitation as for the 2D-SPM (i.e., the characteristically higher Br concentration during the second irrigation [Fig. 8] is not captured), since both models assume a single mobile pore domain. Field sites where flow and transport through the SM domain cannot be neglected require a two-domain concept that allows for flow in both PF and SM domains. Application of MIM-type DPM is limited and clearly results in unrealistic solute mixing volumes; some boundary conditions for solutes cannot be implemented. Moreover, the fraction of the mobile pore region (mobile water content parameter
sm) from an adequate calibration of flow parameters is often too large (questioning the physical meaning of effective parameters; see Haws et al. [2005]) as required for simulating rapid arrival times and when compared with observations of staining patterns. Soil and drainage conditions as reported here (cracked clay and structured macroporous soil) where drainage occurs only during and after storms seem to be widespread for tile-drained field with glacial till or clay soils (Villholth and Jensen, 1998; Villholth et al., 1998; Gärdenäs et al., 2006) at many other sites, including the northern parts of the USA. For other tile-drained field soils, showing "weak" structures (e.g., Infeld site; Köhne et al., 2006) or field soil conditions with gravel at the bottom (e.g., Weiherbach catchment; Zehe and Flühler, 2001), it may be less important. Such soil structure–related effects may explain contrasting interpretations of the results.
Simplifying Hydraulic Assumptions
Furthermore, the simulations still use the drainage branch of the hysteretic water retention curves for simulating water infiltration conditions. The fit between measured and simulated drainage hydrographs (Fig. 3) seems to improve with the duration of relatively wet soil conditions, indicating that for dryer periods, other processes not accounted for in the models became relatively more important (i.e., runoff on crusted or frozen soil surfaces, interception losses, topography and lane effects on local ponding patches, neglected deep seepage rates, or groundwater flow parallel to the tile lines). Independently estimated DPERM parameters (Köhne et al., 2002b) could not be evaluated since they caused numerical problems in conjunction with the boundary conditions assumed here. Also, measurement problems, especially of precipitation during winter (Peck, 1997) may have contributed to some of the smaller discrepancies (Fig. 3). For a better understanding of complex mechanisms involved in preferential flow and transport toward drain tiles, the modeling and experimental challenge is how to determine the upper boundary conditions to represent the infiltration, separated for soil domains or depending on distributed and transient changing soil hydraulic properties during storms or irrigation.
Potassium Bromide Salt Application Effect
Tracer experiments at the Bokhorst site and similar glacial till sites mostly reported that Br effluent concentrations are highest in the first drainage event following solute application and are successively decreasing in subsequent flow events (Gerke and Köhne, 2004; Villholth et al., 1998; Kung et al., 2000b). In this plot scale experiment (Fig. 11
), the Br concentration was largest for the second day of irrigation. Furthermore, there was no dilution-induced concentration drop on the second day but a stronger and faster increase of Br concentrations in tile-drain effluent compared with the first irrigation. The comparison of drainage Br effluent curves (Fig. 11) demonstrates the importance of the second irrigation for initiating the tracer flush toward the drain and that of the following rainfall event for dilution of Br concentrations at the field scale. The measured data indicate a small Br background concentration level, not included in the simulations; thus, the initial concentration drop (Day 97) could not be simulated (Fig. 11, bottom).

View larger version (29K):
[in this window]
[in a new window]
|
FIG. 11. . Measured time series of drain discharge and Br concentration in tile effluent (top) for the irrigation period and the following rain event (Days 97–103), and (bottom) measured compared with simulation results for the irrigation period only (Days 97–99.5). 2D-DPERM, two-dimensional dual-permeability model.
|
|
The main reason for the unexpectedly larger Br effluent concentration peak of the second irrigation is probably related to the application of KBr salt on the surface of the relatively moist soil profile instead of spraying a liquid Br solution. The results for the two scenarios (Fig. 7 and 8) suggest that the salt dissolution process during irrigation forced most of the Br solution to infiltrate into the SM domain at the soil surface. Hence, a short time after the start of irrigation, the kinetically controlled Br transfer from the SM into the PF domain was restricting Br leaching. Since in the topsoil, both domains rapidly approached similarly high moisture contents, the "convective" transfer component was reduced such that mass transfer was then dominated by diffusion, which generally proceeds more slowly but for a longer period of time. On the second day, higher Br effluent concentrations appeared since the SM domain contained already increased levels of Br concentrations in soil water solution at greater depths that could enter the PF domain and be leached more rapidly.
The simulation results (Fig. 7 and 8) demonstrate the considerable effect of the boundary conditions in the 2D-DPERM, which has been rarely discussed previously. Comparing Fig. 7 with Fig. 8 indicates that optimal simulation results would have been obtained when changing from SM application at Day 97 to application in SM and PF domain at Day 98. The physical basis for this change in boundary conditions is the effect of the rate-limited dissolution of the applied KBr salt. The KBr slowly dissolved during the application on Day 97 and initially released only small concentrations to the fast infiltrating bypass water in the PF domain, while presumably being at equilibrium with the slowly infiltrating matrix solution. On Day 98, when all salt was completely dissolved, high concentrations were available to both domains. Given the significant effect on simulation results, the proper treatment of boundary conditions in DPERM approaches deserves more attention. Furthermore, research is necessary to adequately measure and simulate temporal change of solute boundary conditions for application of solid salts (e.g., mineral fertilizers) or of other agrochemicals applied as suspensions with concentrations above solubility (e.g., many pesticides). The uncertainties involved in the soil solute boundary condition and solute mass transfer rate parameters are reflected in the quite variable Nash–Sutcliffe model efficiency coefficients (Table 3).
Two-Dimensional Flow Field Effects
The 2D flow field is strongly affected by surface infiltration and rainfall rates used in the simulations. Temporary ponding at the soil surface during rain or irrigation may occur. Here, for the irrigation of 4.5 h and 3.5 h, respectively, averaged intensities were used for 2D-SPM simulation while real "0.5-h switching" intensities were used in the 2D-DPERM runs. In case of 2D-SPM, the use of the "real" irrigation scheme would have resulted in massive saturation and water balancing problems. By contrast, in the 2D-DPERM simulations, ponding on the surface of the SM domain could be accounted for and redirected to infiltrate into the PF domain (i.e., added to atmospheric BC). However, this representation of the infiltration process, while being more realistic than in previous modeling attempts (Köhne and Gerke, 2005), is still highly simplified compared with natural conditions. Similarly, as it was not possible to account for the local heterogeneity of soil hydraulic properties (except for layering and domains), it was beyond our scope to consider the scattered local patches of ponded water as observed in spots of low topography and within a compacted tractor lane passing through one end of the irrigated partial area without Br application.
Lateral Bromide Plume Displacement
The 2D-SPM and 2D-DPERM simulations could not adequately reproduce the lateral spreading of Br concentrations (Fig. 9). The measured Br plume spread more to the opposite direction than to the drain tile, probably caused by perched water table on a plow pan. The measured water table could not be simulated with the 2D-DPERM using scenarios without parameter calibration. As a result, however, the smearing of concentration fronts due to dilution, mass transfer, and mixing of the solutes in the field was probably larger compared with simulations. Additional "smearing" possibly resulted from effects of spatial heterogeneity of soil properties and temporal changes due to swelling and shrinking, which were not taken into account here. Moreover, the refilled soil material surrounding the tiles could have imposed additional hydraulic resistances that may have restricted drainage. Tensiometer data (e.g., Fig. 5) indicated generally higher water saturation compared with the simulated pressure head profiles. Nevertheless, this first attempt to compare simulations with both the response at the tile-drain outlet as well as with field data of resident Br concentrations is promising since model interpretations do not seem to be inconsistent. The study confirms previous conclusions that model-based analyses of field tracer experiments solely based on outflow data are not sufficient. Drainage curves may be simulated comparably with different model approaches. Since effluent data are integrated over the catchment, both the larger (catchment) scale effects and boundary effects need to be considered, while at the same time, the local effects need to be conserved to regain the residual distribution.
Plot- and Field-Scale Mixing Effect
Although a 3D approach would be better in principle to combine plot scale results with those from other parts of the drain catchments, the advantage of a 3D approach would be limited by our restricted knowledge on boundary conditions and spatial distribution of soil hydraulic properties at this site and for the time being appeared beyond the scope of this work. Note that we did not intend to optimize the match between simulated and measured values by adjusting parameters. Of course, a better resemblance would be possible especially when including plow pan and other effects but comparisons would become more complex.
Effects of the complex 3D subsurface flow field (e.g., Troch et al., 2003) on Br leaching are difficult to evaluate in detail. For example, the rainfall event following the irrigation did not stimulate Br leaching, although tile flow rates strongly increased (Fig. 11). As a result of the irrigation, soil moisture and water table at the plot were already increased compared with the rest of the field, and Br has been spread in the soil before the rain. Furthermore, rainfall intensities were lower than local irrigation rates, which reduced the bypass component of preferential Br transport. Thus, drainage from the plot started faster than before and caused dilution of Br concentrations in drain effluent. At the irrigated plots, the moisture regime was then quite different from that of the surrounding nonirrigated areas. Consequently, the drainage contributions from irrigated and nonirrigated areas to the composite tile outflows differed in a highly dynamic way for a few days after irrigation, after which this "partial-area irrigation" effect faded away. In the simulations we tried to take into account water contributions from within vs. outside the plot by separately simulating 2D representations of plots with and without irrigation and applying a plot- and field-scale weighing procedure (see Appendix). Still, a dilution effects from drainage contribution of the irrigated plot of 120 m2 was more exaggerated in the simulations where the Br concentrations started to rise at Day 100, then sharply dropped and later increased again for 2D-DPERM (Fig. 8). Here, the resolution of measured effluent concentrations was probably not sufficient to clearly reveal these plot scale effects during Day 100. Interestingly, the drop in measured Br effluent concentrations during Day 97 (Fig. 8, top) could be interpreted similarly. Here, the background Br effluent concentration decreased probably because of dilution by irrigation-induced drain water that arrived shortly before the preferentially transported applied Br appeared at the tile outlet. However, such interpretations based on few observation and small concentration differences remain speculative without further experimental verification.
 |
Conclusions
|
|---|
A strongly improved description of preferential flow toward tile drain is obtained when using the 2D-DPERM compared with the 2D-SPM and previously published 1D descriptions. The same hydraulic parameters were used here as in a 1D-DPERM study for another tracer experiment at the same Bokhorst field site. Although parameters were not calibrated, and a number of effects including soil spatial heterogeneity were not considered, this 2D-DPERM captured most of the basic processes at three scales (structure, pedon, and plot), as well as the dimensionality of the problem: (i) profile-scale bypass-flow and physical nonequilibrium, (ii) 2D lateral flow toward a drain, and (iii) 3D field water balance although in simplistic form (plot-scale and field-scale mixing effects). The superiority of 2D-DPERM over 2D-SPM was further supported by the consistency in which simulation results compared better with both the response at the tile drain outlet as well as the resident Br concentrations. A comparison with MIM results was not possible using the same parameters and scenarios that assume SM domain involvement in macroscopic flow and transport processes.
The study confirms previous conclusions that model-based analyses of field tracer experiments solely based on outflow data are not sufficient; local effects need to be conserved to regain the residual distributions. Spatially distributed fluxes in the unsaturated and saturated soil are needed for an improved understanding of complex mechanisms involved in preferential flow and transport toward subsurface drains. Results of the simulated irrigation events indicate that the soil surface solute boundary conditions in the 2D-DPERM had significant effects on tile effluent leaching patterns. The 2D-DPERM simulations suggest that most of the applied KBr salt that dissolved during the irrigation first entered the SM domain before being transferred into the PF domain. Later, probably both domains contributed to Br infiltration. Given the significant effect on simulation results, the proper treatment of boundary conditions in DPERM approaches deserves more attention. The characteristic feature of a higher Br concentration in tile effluent during the second irrigation was captured here with less-restricted solute mass transfer as assumed in previous 1D-DPERM simulations. Due to the half-hour alternating irrigation scheme, the SM domain could have stored relatively more Br while more water infiltrated the PF domain during irrigation. Increased soil moisture at the irrigated plot has provoked more rapid drainage during following rainfall event than in the other areas, as well as more rapid Br leaching. Although simplified, the plot- and field-scale weighing procedure accounted for increased plot scale soil moisture effects on total flow and also for dilution effects for Br effluent calculation (i.e., 60 m2 vs. 5000 m2).
In contrast to previous 1D-DPERM analyses where mass transfer effects were found most sensitive, the 2D-DPERM analysis suggests that boundary conditions at the soil surface, near the water table, and of the field-scale mixing are significantly affecting leaching patterns in addition to physical nonequilibrium effects. Further improvements of DPERM analyses require developments in numerical modeling and field experimentation, including special attention toward spatially distributed and temporally dynamic soil structure and soil surface conditions.
 |
Appendix: Weighing Procedures for Approximating Field-Scale Tile Drainage and Leaching from Plot Scale Simulations
|
|---|
Weighing from DPERM concept yields the specific drain discharge (water fluxes in centimeters per day), the specific drain discharge from the plot (Eq. [1c]) is
 | [A1] |
Weighing from the 2D vertical cross-section concept, drainage flux, qdrain, per unit length of the seepage face (i.e., Ldrain = 4 cm) is converted into specific discharge, qsurf (per unit length of the cross-section surface (i.e., Lsurf = 590 cm) by
 | [A2] |
Weighing of surface irrigated area in relation to the tile drain catchment area yields
 | [A3] |
The total catchment drain discharge, qdrain, is obtained as
 | [A4] |
Br concentration of drained solutes related to application area (mg cm–3) for 2D-SPM is calculated as
 | [A5] |
and for 2D-DPERM as
 | [A6] |
The effluent concentration in mg L–1 related to total catchment area, ctotal, is
 | [A7] |
 |
ACKNOWLEDGMENTS
|
|---|
We thank the Ministry of Education of the Czech Republic for financial support under grant MSM 6840770002, the German Federal Ministry of Consumer Protection, Food and Agriculture under grant 05/04 (BMELV), and Ministry for Rural Development, Environment and Consumer Protection (MLUV) of the state of Brandenburg. Additional support was provided through the bilateral German–Czech research cooperation (KONTAKT 05/04). The comments of Dr. S. Köhne (University of Rostock) and two anonymous reviewers are gratefully acknowledged.
 |
REFERENCES
|
|---|
- Arabi, M., J.S. Stillman, and R.S. Govindaraju. 2006. A process-based transfer function approach to model tile-drain hydrographs. Hydrol. Proc. 20:3105–3117.[CrossRef]
- Bear, J. 1972. Dynamics of fluids in porous media. American Elsevier, New York.
- Boivin, A., J.
im
nek, M. Schiavon, and M.Th. van Genuchten. 2006. Comparison of pesticide transport processes in three tile-drained field soils using HYDRUS-2D. Vadose Zone J. 5:838–849.[Abstract/Free Full Text] - Bosch, D., K. King, R. Yoder, J. Parsons, W. Rawls, L. Bergstrom, G. Brown, R. Malone, A. Ward, A. Shirmohammadi, and G. Vellidis (ed.). 2000. Preferential flow, water movement and chemical transport in the environment. Proc. ASAE 2nd Int. Sym., Honolulu, HI. 3–5 Jan. 2001. ASAE, St. Joseph, MI.
- Bronswijk, J.J.B., W. Hamminga, and K. Oostindie. 1995. Field-scale solute transport in a heavy clay soil. Water Resour. Res. 31:517–526.[CrossRef]
- Czapar, G.F., R.S. Kanwar, and R.S. Fawcett. 1994. Herbicide and tracer movement to field drainage tiles under simulated rainfall conditions. Soil Tillage Res. 30:19–32.[CrossRef]
- Elliott, J.A., A.J. Cessna, W. Nicholaichuk, and L.C. Tollefson. 2000. Leaching rates and preferential flow of selected herbicides through tilled and untilled soil. J. Environ. Qual. 29:1650–1656.[Web of Science]
- FAO. 1998. World reference base for soil resources. World Soil Resources Report No. 84. ISSS, ISRIC, FAO, Rome, Italy.
- Fox, G.A., G.J. Sabbagh, W.L. Chen, and M.H. Russell. 2006. Uncalibrated modelling of conservative tracer and pesticide leaching to groundwater: Comparison of potential Tier II exposure assessment models. Pest Manage. Sci. 62:537–550.[CrossRef]
- Feddes, R.A., P.J. Kowalik, and H. Zaradny. 1978. Simulation of field water use and crop yield. PUDOC, Wageningen, The Netherlands.
- Flury, M. 1996. Experimental evidence of transport of pesticides through field soils: A review. J. Environ. Qual. 25:24–45.
- Fortin, J., E. Gagnon-Bertrand, L. Vézina, and M. Rompré. 2002. Preferential bromide and pesticide movement to tile drains under different cropping practices. J. Environ. Qual. 31:1940–1952.[Abstract/Free Full Text]
- Gärdenäs, A.I., J.
im
nek, N. Jarvis, and M.Th. van Genuchten. 2006. Two-dimensional modelling of preferential water flow and pesticide transport from a tile-drained field. J. Hydrol. 329:647–660.[CrossRef] - Gaur, A., R. Horton, D.B. Jaynes, and J.L. Baker. 2006. Measured and predicted solute transport in a tile drained field. Soil Sci. Soc. Am. J. 70:872–881.[Abstract/Free Full Text]
- Gerke, H.H. 2006. Preferential flow descriptions for structured soils. J. Plant Nutr. Soil Sci. 169:382–400.[CrossRef]
- Gerke, H.H., and J.M. Köhne. 2002. Estimating hydraulic properties of soil aggregate skins from sorptivity and water retention. Soil Sci. Soc. Am. J. 66:26–36.[Abstract/Free Full Text]
- Gerke, H.H., and J.M. Köhne. 2004. Dual-permeability modeling of preferential bromide leaching from a tile drained glacial till agricultural field. J. Hydrol. 289:239–257.[CrossRef]
- Gerke, H.H., and M.Th. van Genuchten. 1993a. A dual-porosity model for simulating the preferential movement of water and solutes in structured porous media. Water Resour. Res. 29:305–319.[CrossRef]
- Gerke, H.H., and M.Th. van Genuchten. 1993b. Evaluation of a first-order water transfer term for variably saturated dual-porosity flow models. Water Resour. Res. 29:1225–1238.[CrossRef]
- Gerke, H.H., and M.Th. van Genuchten. 1996. Macroscopic representation of structural geometry for simulating water and solute movement in dual-porosity media. Adv. Water Resour. 19:343–357.[CrossRef]
- Haude, W. 1955. Zur Bestimmung der Verdunstung auf möglichst einfache Weise. (In German.) Mitt. d. DWD 11:1–24.
- Haws, N.W., W.P. Ball, and E.J. Bouwer. 2006. Modeling and interpreting bioavailability of organic contaminant mixtures in subsurface environments. J. Contam. Hydrol. 82:255–292.[CrossRef][Web of Science][Medline]
- Haws, N.W., B.S. Das, and P.S.C. Rao. 2004. Dual-domain solute transfer and transport processes: Evaluation in batch and transport experiments. J. Contam. Hydrol. 75:257–280.[CrossRef][Web of Science][Medline]
- Haws, N.W., P.S.C. Rao, J.
im
nek, and I.C. Poyer. 2005. Single-porosity and dual-porosity modeling of water flow and solute transport in subsurface-drained fields using effective field-scale parameters. J. Hydrol. 313:257–273.[CrossRef] - Jarvis, N.J. 1994. The MACRO model, version 3.1: Technical description and sample simulations. Rep. and diss. 19. Dep. of Soil Science, Swedish Univ. of Agricultural Sciences, Uppsala.
- Jaynes, D.B., S.I. Ahmed, K.-J.S. Kung, and R.S. Kanwar. 2001. Temporal dynamics of preferential flow to a subsurface drain. Soil Sci. Soc. Am. J. 65:1368–1376.[Abstract/Free Full Text]
- Jaynes, D.B., K.-J.S. Kung, S.I. Ahmed, and R.S. Kanwar. 2000. Temporal dynamics of preferential flow to a field tile. p. 89–92. In D. Bosch, K. King, R. Yoder, J. Parsons, W. Rawls, L. Bergstrom, G. Brown, R. Malone, A. Ward, A. Shirmohammadi, and G. Vellidis (ed.) Preferential flow, water movement, and chemical transport in the environment. Proc. ASAE 2nd Int. Symp., Honolulu, HI. 3–5 Jan. 2001. ASAE, St. Joseph, MI.
- Kamra, S.K., J. Michaelsen, W. Wichtmann, and P. Widmoser. 1999. Preferential solute movement along the interface of soil horizons. Water Sci. Technol. 40:61–68.[Web of Science]
- Kasteel, R., M. Burkhardt, S. Giesa, and H. Vereecken. 2005. Characterization of field tracer transport using high-resolution images. Vadose Zone J. 4:101–111.[Abstract/Free Full Text]
- Kladivko, E.J., L.C. Brown, and J.L. Baker. 2001. Pesticide transport to subsurface tile drains in humid regions of North America. Crit. Rev. Environ. Sci. Technol. 31:1–61.[CrossRef]
- Kohler, A., K.C. Abbaspour, M. Fritsch, and R. Schulin. 2003. Using simple bucket models to analyze solute export to subsurface drains by preferential flow. Vadose Zone J. 2:68–75.[Abstract/Free Full Text]
- Kohler, A., K.C. Abbaspour, M. Fritsch, M.Th. van Genuchten, and R. Schulin. 2001. Simulating unsaturated flow and transport in a macroporous soil to tile drains subject to an entrance head: Model development and preliminary evaluation. J. Hydrol. 254:67–81.[CrossRef]
- Köhne, J.M. 1999. Analyse präferentiellen Wasserflusses und Stofftransports in strukturierten Böden mit Hilfe eines Dual-Porositätsmodells. Ph.D. diss. (In German, with English summary.) In H. Roweck and B. Lennartz (ed.) Schriftenreihe Inst. Water Management and Landscape Ecology 30. Christian-Albrechts-Univ., Kiel, Germany.
- Köhne, J.M., and H.H. Gerke. 2005. Spatial and temporal dynamics of preferential tracer movement towards a tile drain. Vadose Zone J. 4:79–88.[Abstract/Free Full Text]
- Köhne, J.M., H.H. Gerke, and S. Köhne. 2002a. Effective diffusion coefficients of soil aggregates with surface skins. Soil Sci. Soc. Am. J. 66:1430–1438.[Abstract/Free Full Text]
- Köhne, J.M., S. Köhne, and H.H. Gerke. 2002b. Estimating the hydraulic functions of dual-permeability models from bulk soil data. Water Resour. Res. 38:1121, doi:10.1029/2001WR000492.
- Köhne, S., B. Lennartz, J.M. Köhne, and J.
im
nek. 2006. Bromide transport at a tile-drained field site: Experiment, and one- and two-dimensional equilibrium and non-equilibrium numerical modeling. J. Hydrol. 321:390–408.[CrossRef] - Kronvang, B., H.L. Strom, C.C. Hoffmann, A. Laubel, and N. Friberg. 2004. Subsurface tile drainage loss of modern pesticides: Field experiment results. Water Sci. Technol. 49:139–147.[Web of Science]
- Kung, K.-J.S., E.J. Kladivko, T.J. Gish, T.S. Steenhuis, G. Bubenzer, and C.S. Helling. 2000a. Quantifying preferential flow by breakthrough of sequentially applied tracers: Silt loam soil. Soil Sci. Soc. Am. J. 64:1296–1304.[Abstract/Free Full Text]
- Kung, K.-J.S., T.S. Steenhuis, E.J. Kladivko, T.J. Gish, G. Bubenzer, and C.S. Helling. 2000b. Impact of preferential flow on the transport of adsorbing and non-adsorbing tracers. Soil Sci. Soc. Am. J. 64:1290–1296.[Abstract/Free Full Text]
- Larsbo, M., S. Roulier, F. Semo, R. Kasteel, and N. Jarvis. 2005. An improved dual-permeability model of water flow and solute transport in the vadose zone. Vadose Zone J. 4:398–406.[Abstract/Free Full Text]
- Larsson, M.H., and N.J. Jarvis. 1999. Evaluation of a dual-porosity model to predict field-scale solute transport in a macroporous soil. J. Hydrol. 215:153–171.[CrossRef]
- Lennartz, B., J. Michaelsen, W. Wichtmann, and P. Widmoser. 1999. Time variance analysis of preferential solute movement at a tile-drained field site. Soil Sci. Soc. Am. J. 63:39–47.[Abstract/Free Full Text]
- Millington, R.J., and J.P. Quirk. 1961. Permeability of porous soils. Trans. Faraday Soc. 57:1200–1207.[CrossRef]
- Mohanty, B.P., R.S. Bowman, J.M.H. Hendrickx, J.
im
nek, and M.Th. van Genuchten. 1998. Preferential transport of nitrate to a tile drain in an intermittent-flood-irrigated field: Model development and experimental evaluation. Water Resour. Res. 34:1061–1076.[CrossRef] - Nash, J.E., and J.V. Sutcliffe. 1970. River flow forecasting through conceptual models part I: A discussion of principles. J. Hydrol. 10:282–290.[CrossRef]
- Peck, E.L. 1997. Quality of hydrometeorological data in cold regions. J. Am. Water Resour. Assoc. 33:125–134.[CrossRef]
- Ray, C., T. Vogel, and J. Dusek. 2004. Modeling depth-variant and domain-specific sorption and biodegradation in dual-permeability media. J. Contam. Hydrol. 70:63–87.[CrossRef][Web of Science][Medline]
- Richard, T.L., and T.S. Steenhuis. 1988. Tile drain sampling of preferential flow on a field scale. J. Contam. Hydrol. 3:307–325.[CrossRef]
- Schelde, K., L.W. de Jonge, C. Kjaergaard, M. Laegdsmand, and G.H. Rubaek. 2006. Effects of manure application and plowing on transport of colloids and phosphorus to tile drains. Vadose Zone J. 5:445–458.[Abstract/Free Full Text]
im
nek, J., N.J. Jarvis, M.Th. van Genuchten, and A. Gärdenäs. 2003. Review and comparison of models for describing nonequilibrium and preferential flow and transport in the vadose zone. J. Hydrol. 272:14–35.[CrossRef]- Stamm, C., H. Fluhler, R. Gachter, J. Leuenberger, and H. Wunderli. 1998. Preferential transport of phosphorus in drained grassland soils. J. Environ. Qual. 27:515–522.[Web of Science]
- Stamm, C., R. Sermet, J. Leuenberger, H. Wunderli, H. Wydler, H. Flühler, and M. Gehre. 2002. Multiple tracing of fast solute transport in a drained grassland soil. Geoderma 109:245–268.[CrossRef][Web of Science]
- Stone, W.W., and J.T. Wilson. 2006. Preferential flow estimates to an agricultural tile drain with implications for glyphosate transport. J. Environ. Qual. 35:1825–1835.[Abstract/Free Full Text]
- Troch, P.A., C. Paniconi, and E.E. van Loon. 2003. Hillslope-storage Boussinesq model for subsurface flow and variable source areas along complex hillslopes: 1. Formulation and characteristic response. Water Resour. Res. 39(11):1316, doi:10.1029/2002WR001728.[CrossRef]
- Vanderborght, J., and H. Vereecken. 2007. Review of dispersivities for transport modeling in soils. Vadose Zone J. 6:29–52.[Abstract/Free Full Text]
- van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.[Web of Science]
- van Genuchten, M.Th., and F.N. Dalton. 1986. Models for simulating salt movement in aggregated field soils. Geoderma 38:165–183.[CrossRef][Web of Science]
- Vereecken, H. 2005. Mobility and leaching of glyphosate: A review. Pest Manage. Sci. 61:1139–1151.[CrossRef]
- Villholth, K.G., and K.H. Jensen. 1998. Flow and transport processes in a macroporous subsurface-drained glacial till soil: II. Model analysis. J. Hydrol. 207:121–135.
- Villholth, K.G., K.H. Jensen, and J. Fredericia. 1998. Flow and transport processes in a macroporous subsurface-drained glacial till soil: I. Field investigations. J. Hydrol. 207:98–120.
- Vogel, T. 1987. SWMS II: Numerical model of two-dimensional flow in a variably saturated porous medium. Research Report No. 87, Dep. of Hydraulics and Catchment Hydrology, Agriculture Univ., Wageningen, The Netherlands.
- Vogel, T., and M. Císlerova. 1988. On the reliability of unsaturated hydraulic conductivity calculated from the moisture retention curve. Transp. Porous Media 3:1–15.[CrossRef]
- Vogel, T., H.H. Gerke, R. Zhang, and M.Th. van Genuchten. 2000. Modeling flow and transport in a two-dimensional dual-permeability system with spatially variable hydraulic properties. J. Hydrol. 238:78–89.[CrossRef]
- Vogel, T., M.Th. van Genuchten, and M. Cislerova. 2001. Effect of the shape of soil hydraulic functions near saturation on variably-saturated flow predictions. Adv. Water Resour. 24:133–144.[CrossRef]
- Wichtmann, W. 1994. Stoffeintrag aus landwirtschaftlichen Dränflächen in Fließgewässer (Messung und Simulation). Ph.D. diss. (In German, with English summary.) Institut für Wasserwirtschaft und Landschaftsökologie der Christian-Albrechts-Universität Kiel, Kiel, Germany.
- Wichtmann, W., B. Lennartz, and P. Widmoser. 1998. Bromidverlagerung an zwei gedränten Standorten in Schleswig-Holstein. (In German.). J. Plant Nutr. Soil Sci. 161:121–128.
- Zadoks, J.C., T.T. Chang, and C.F. Konzak. 1974. A decimal code for the growth stages of cereals. Weed Res. 14:415–421.[CrossRef]
- Zehe, E., and H. Flühler. 2001. Preferential transport of isoproturon at a plot scale and a field-scale tile-drained site. J. Hydrol. 247:100–115.[CrossRef]
This article has been cited by other articles:

|
 |

|
 |
 
J. Dusek, H. H. Gerke, and T. Vogel
Surface Boundary Conditions in Two-Dimensional Dual-Permeability Modeling of Tile Drain Bromide Leaching
Vadose Zone J.,
November 1, 2008;
7(4):
1287 - 1301.
[Abstract]
[Full Text]
[PDF]
|
 |
|