Quantifying the impact of preferential flow on solute transport to tile drains in a sandy field soil

HYDROL 3683

Journal of Hydrology 215 (1999) 116–134

Quantifying the impact of preferential flow on solute transport to tile drains in a sandy field soil M.H. Larsson a,*, N.J. Jarvis a, G. Torstensson a, R. Kasteel b a

Department of Soil Sciences, Swedish University of Agricultural Sciences, Box 7014, 75007 Uppsala, Sweden b Department of Environmental Physics, University of Heidelberg, D-69120, Heidelberg, Germany Received 15 December 1997; accepted 30 June 1998

Abstract The objective of this study was to quantify the impact of preferential finger flow on solute leaching to tile drains and shallow groundwater in a water repellent sandy soil. Measurements of site hydrology and bromide movement were made during a sixmonth period following application in autumn 1994 to a 900 m 2 plot. A dye tracing study was also conducted and this confirmed the presence of finger flow in the water repellent topsoil. Water flow and bromide transport was simulated using a mobileimmobile concept, in which a fraction of the soil volume is assumed not to participate in transport. With calibration of ‘difficult’ parameters (e.g. the fractional mobile volume), and accounting for significant water inflows to the plot from the surroundings, there was a good agreement between simulated and observed hydrology and solute transport patterns at the site. Mass balance calculations showed that ca. 46% of the bromide application was lost to the tile drains in the winter period following application, while 16% was lost in shallow lateral groundwater flow. Comparisons between one- and two-domain simulations showed that preferential flow reduced soil water contents and increased drainflow in the early autumn, leading to a more rapid bromide transport to depth. Nevertheless, the impact of preferential finger flow on leaching was minor for the non-reactive tracer bromide in this experiment, presumably because it was applied to relatively wet soil. However, scenario simulations for autumn application of a short half-life pesticide showed more significant effects, with leaching increased by ca. 80% because of preferential flow (from 1.2% to 2.2% of the applied amount). 䉷 1999 Elsevier Science B.V. All rights reserved. Keywords: Field-scale; Water repellency; Dye tracing; Finger flow; Solute transport; Mobile–immobile model

1. Introduction Preferential flowpaths may allow rapid movement of solutes from the soil surface to the groundwater, by-passing a significant volume of the unsaturated zone. The interaction between the soil matrix and potential pollutants such as heavy metals, pesticides, nitrate and phosphate through adsorption, immobilization, and degradation processes is diminished, * Corresponding author. Fax: ⫹ 46-18-673430; e-mail: martin. [email protected]

thereby increasing the risk of leaching of these solutes to the groundwater or to surface waters via subsurface drains. Indeed, contaminants have often been found to move to greater depths at faster rates than would be predicted by transport theory based on continuum approaches applied on the basis of area-averaged water potentials and solute concentrations (see reviews by Jury and Flu¨hler, 1992; Jarvis, 1998). These preferred flow paths may be of different types. Beven and Germann (1982) defined macropores as structural pores (e.g. biotic pores formed by earthworms and plant roots, cracks produced by

0022-1694/99/$ - see front matter 䉷 1999 Elsevier Science B.V. All rights reserved. PII: S0022-169 4(98)00265-0

M.H. Larsson et al. / Journal of Hydrology 215 (1999) 116–134

shrinkage of clay soils, chemical weathering of bedrock material, freeze/thaw cycles and soil cultivation) that permit rapid non-equilibrium channel flow, irrespective of size. Preferential flow can also occur in the matrix of sandy soils with no obvious structural channels. The underlying mechanisms may vary and include air pressure build-up ahead of the wetting front, water repellency (e.g. Bond, 1964; Dekker and Jungerius, 1990; Hendrickx et al., 1993) and spatial variation in soil texture, for example, at horizon boundaries where a finer soil overlies a coarse, dry sand layer (e.g. Hill and Parlange, 1972; Glass et al., 1990), or where inclined coarse sand lenses interbedded in a finer-textured soil matrix induce lateral redistribution of water to form concentrated ‘column’ flows (Kung, 1990). Development of mechanistic models which can describe preferential flow in sandy soils under field conditions is still at a very early stage, although some attempts have recently been made (e.g. Hendrickx and Yao, 1996; de Rooij, 1996; Nieber, 1996). Much of the theoretical work on unstable flow has concentrated on understanding the processes as they occur in ‘ideal’ porous media in the laboratory (Glass et al., 1989). It is still unclear to what extent this experience and knowledge can be transferred to the field situation, so that the conditions under which finger flow can be expected in the field cannot be easily predicted ‘a priori’ and nor can its consequences (Glass and Nicholl, 1996). Thus, although model parameters must be determined from measurements by calibration or ‘inverse’ modelling, functional approaches such as the mobile/immobile water convection-dispersion equation (MIMCDE) may still prove useful tools to analyse preferential flow and transport in the matrix of coarsertextured field soils. Analytical solutions of the MIMCDE have been developed for both non-reactive and sorbing and degrading solutes (e.g.van Genuchten and Wierenga, 1976; van Genuchten and Wagenet, 1989; Gamerdinger et al., 1990). These analytical solutions require idealised initial and boundary conditions and assume steady state water flow, and are therefore not well-suited to field situations characterized by timevarying soil water contents. For numerical solutions which deal with transient conditions, one unresolved question is how best to define the fractional mobile water content in relation to the time-varying total water content (see discussion in Zu¨rmuhl and Durner,

117

1996). There may be no single best answer to this question: the most appropriate solution may well depend on the particular nature of the flow and transport processes in the soil under study. van Dam et al. (1990; 1996; 1997) presented a simple variant of the MIM-CDE concept which may be particularly appropriate to water repellent soils where a given fraction of the soil matrix remains effectively unwetted during the infiltration process. In this paper, we used a modified version of the MIM-CDE model introduced by van Dam et al. (1990) to quantify the effects of preferential water flow on non-reactive solute leaching (Br ⫺) measured in tile-drained plots in a water repellent sandy soil. Quantifying the impact of preferential flow on nonreactive solute leaching depends entirely on the time frame for comparison, because all the dose must sooner or later leach from the profile. Therefore, scenario simulations were also run using the calibrated model to quantify the effects of preferential flow on leaching of a hypothetical reactive solute representing a short half-life, mobile pesticide.

2. Materials and methods 2.1. Site details The experiments were carried out at the Mellby field site (56⬚29 0 N) in south-west Sweden, which consists of tile-drained plots, each 30 by 30 m in size. Plot 39 was used for the tracer experiment described in this paper, while the neighbouring plot 40 was used for more destructive complementary studies, including a dye tracing experiment, and soil sampling for hydraulic properties. The primary drainage system at the site consists of perforated plastic pipes at a spacing of 6.75 m and at a depth of ca. 90 cm. Each plot is also surrounded by a single discard drainage pipe at the same depth as the primary drainage system. This secondary drain, which was not monitored, is intended to minimize lateral flows between plots. The plots have been in continuous spring cereal cultivation since 1993. Meteorological data were obtained on site with a standard automatic weather station recording precipitation, air temperature, wind speed, vapour pressure and global radiation on an hourly basis.

118

M.H. Larsson et al. / Journal of Hydrology 215 (1999) 116–134

Table 1 Particle size distribution and organic matter contents in plot 39 at Mellby (n ˆ 8, standard deviations in parentheses) Depth (cm)

Sand (%) 60–2000 mm

Fine sand (%) 60–200 mm

Silt (%) 2–60 mm

Clay (%) ⬍ 2 mm

O.M. (%)

0–10 cm 40–50 cm 70–80 cm

84.6 (1.4) 91.4 (4.3) 87.4 (4.2)

36.1 (2.2) 38.9 (3.9) 50.8 (4.8)

8.0 (0.6) 5.1 (3.2) 10.6 (4.2)

7.4 (0.9) 3.5 (2.0) 2.0 (0.6)

4.3 (0.2) 1.2 (0.7) a 0.7 (0.2)

a

Excluding one sample with 6.2% O.M., reflecting infilled topsoil material above a drain trench.

The soil profile at the experimental site consists of three different horizons (Bergstro¨m et al., 1994). The topsoil is a loamy sand in texture. It has a brownish black colour, a large organic matter content and a weak granular structure. The topsoil thickness at Mellby can vary considerably over short distances. A total of 81 corings in plot 39 gave an average topsoil thickness of 32 cm, with a standard deviation of 4 cm, a minimum value of 20 cm and a maximum of 43 cm. The boundary between the topsoil and the subsoil is very sharp. The subsoil B horizon is a yellowish brown structureless sand of low organic matter content. There is also a clay-textured C horizon at Mellby, generally deeper than 1 m in the profile. Unlike the topsoil and the subsoil, nothing is known about the physical properties of this clay layer.

Fig. 1. Soil water retention data (Wind evaporation method).

Fig. 2. Unsaturated hydraulic conductivity, (a): 0–20 cm depth and (b): 40–60cm depth.

M.H. Larsson et al. / Journal of Hydrology 215 (1999) 116–134

119

2.2. Soil physical and hydraulic properties

2.3. Water repellency

Table 1 shows the particle size distributions and organic matter contents measured on eight replicate cores at three depths (0–10 cm, 40–50 cm and 70– 80 cm). Table 1 shows that the soil texture at Mellby appears to be rather homogeneous, with very small coefficients of variation (CV) (generally ⬍5%) for the dominating particle size classes. The subsoil at 70–80 cm depth is dominated by the finer sand fraction and has a larger silt content, whereas the upper subsoil contains more medium and coarse sand. Soil hydraulic properties were measured in the topsoil (0–20 cm depth) and in the subsoil at 40– 60 cm depth. In contrast, no measurements of the hydraulic properties in the deeper subsoil are available. The clay layer is completely saturated most of the time, and is presumed to act as an impermeable lower boundary to the system. In the model simulations it is treated as such, with the base of the profile set at 110 cm depth. Near-saturated hydraulic conductivities were measured by tension infiltrometers in the topsoil and the subsoil at 40 cm depth (Jarvis and Messing, 1995, and unpublished data). In an intermediate saturation range, hydraulic conductivity was measured on two replicate 20 cm diameter cores using the steady-state drip infiltration method. In the dry range, soil water retention curves and unsaturated hydraulic conductivity were measured simultaneously using the Wind evaporation method (Wind, 1966; Stolte et al., 1992). The measurements of soil hydraulic properties, made using the three different methods, are shown in Figs. 1 and 2a,b. One striking feature of the results, particularly from the Wind evaporation method, is the rather large variability in unsaturated hydraulic conductivity (up to 1 order of magnitude) observed in what is apparently a texturally-homogeneous soil (Table 1). The link between variability in soil texture and hydraulic properties and its effects on the pattern of solute transport at Mellby is the subject of ongoing research, but is beyond the scope of the present paper. Here, we investigate the ability of a ‘functional’ deterministic model to describe observed plot-scale water flow and bromide transport using effective soil parameters.

During the dye experiment reported in this paper, water droplets were pipetted onto the soil surface and remained there for several minutes before infiltrating. The organic matter content in the topsoil is quite high (ca. 4%), and although the clay content is not negligible (ca. 7%), soils with clay contents as high as 20% have been reported to be water repellent (McGhie and Posner, 1980). These are all indications that the soil at Mellby is at least potentially water repellent, so that additional laboratory experiments were carried out using the water drop penetration time (WDPT) test to quantify the degree of hydrophobicity. In the WDPT-test, a water droplet is placed on a smoothed air-dry surface of a soil sample and the time to penetrate the soil is measured. Theoretically, a soil with a contact angle equal to, or greater than, 90⬚ should support a water drop on its surface until it evaporates. If the angle is less than 90 ⬚, capillary forces will pull water into the soil. In general, a soil is considered at least initially water repellent if the water drop penetration time exceeds 5 s (van Dam et al., 1990; Dekker and Jungerius, 1990; Hendrickx et al., 1993). The 5 s period was selected for convenience and does not have any specific physical meaning (Richardson, 1984). Disturbed bulk soil samples taken from zero to 20 cm and from 50 to 70 cm depth were first airdried for two weeks. Droplets of distilled water from a pipette were then placed on smooth air-dry soil surfaces and the penetration time was determined. Using this experimental protocol, only a slight water repellency was found for the topsoil, with a mean WDPT value from 84 measurements of 3 s and a maximum value of 11 s. No water repellency was found for the subsoil (WDPT values were always zero). Dekker et al. (1999, this issue) sampled the topsoil at Mellby more extensively and reported somewhat longer WDPT for soil oven-dried at 65⬚C, and a strong depth dependence of WDPT. They found that ca. 95% of samples at 0–5 cm depth had WDPT between 5 and 60 s, with the remainder showing WDPT ⬍ 5 s. Fewer than 30% of samples at 10– 25 cm depth had WDPT between 5 and 60 s, while all samples at 30–35 cm depth had WDPT less than 5 s. Based on these measurements, the upper topsoil at Mellby would be classified as potentially water repellent.

120

M.H. Larsson et al. / Journal of Hydrology 215 (1999) 116–134

Table 2 Soil sampling for volumetric water content and bromide contents Sampling occasion

Number of soil cores

Accumulated precipitation (mm)

Date

Days after application

1 2 Application 3 4 5

16 16

0.0 62.8 62.8 90.9 184.2 426.1

1994-09-22 1994-10-09 1994-10-11 1994-11-08 1994-12-13 1995-03-09

⫺19 ⫺2 0 28 63 149

15 10 14

2.4. Field measurements A mixture of approximately 290 l water and 10 kg of KBr was applied on plots 39 and 40 with a tractor mounted sprayer on 11th October 1994. To measure the mean amount and the spatial variability of the applied solute, 25 aluminium trays, each with an area of 54.7 cm 2 were placed on the surface of plot 39 before spraying. After spraying, the trays were rinsed with a fixed known volume of distilled water and the water was stored at ⫹ 4⬚C until analysed. Br ⫺ analysis showed a mean application of 3.24 g m ⫺2, with a CV of 19% and minimum and maximum individual values of 2.24 and 5.13 g m ⫺2 respectively. Randomly distributed undisturbed soil cylinders were taken on seven occasions during a nearly one year experimental period following application using the motor-driven cylindrical core (1 m long, 9.6 cm internal diameter) described by Hendrickx et al. (1991). However, this paper only presents the results from two samplings prior to application (the first is used as the initial condition in simulations) and for the first three samplings following application (until 9th March 1995, see Table 2). This is because most of the Br ⫺ pulse had already leached from the profile by the end of the winter period. Soil samples were taken from within each cylinder at 10 cm intervals to 90 cm depth using small cores, 203 cm 3 in volume (7.2 cm diameter, 5.0 cm height). Thus, there was no risk of sample cross-contamination, as the outermost layer of the soil within the cylinder was excluded from the samples removed with the small cores. All samples were first weighed and then frozen within 12 h after sampling and stored frozen until analysed. Soil solution samples were extracted from 200 g samples of field-moist soil with 200 ml of deionized water, after shaking for more than 12 h.

Both the soil water extracts and water samples from the drainage pipes were analysed using a Dionex 2000 ion chromatograph. The detection limit for Br ⫺ in water samples was 0.1 mg l ⫺1, and was estimated to be 0.2 g m ⫺3 for Br ⫺ in soil. Gravimetric water content was determined by oven-drying ca. 100 g of soil at 105⬚C for 24 h. Bulk density and volumetric water content were calculated from the known weight of the fresh sample and the core volume. Drainage flow from the plots was measured continuously using a calibrated tipping bucket system, with the number and timing of tips recorded by a data logger. The resolution of the measured drainflow was 0.008 mm. Flow-proportional drainage water samples were taken by an automatic flow sampler and stored in plastic bottles for later analysis. One portion of water was withdrawn for each 0.155 mm of drainflow (as recorded by the tipping buckets), and each bottle contained ten such portions, so that each analysed sample represents 1.55 mm of drainflow. It became apparent from a comparison of drainflow with measured precipitation inputs that considerable amounts of water must have flowed into plot 39 either laterally from the surrounding plots or as artesian groundwater from deeper layers. Indeed, for much of the winter period 1994/5, measured drainage from plot 39 substantially exceeded precipitation. At the same time, drainflow was negligible from several of the surrounding plots, including plot 40. This failure of the drainage systems was presumably caused by clogging of the pipes by iron oxides (the Mellby subsoil is rich in iron oxides), and is thought to be the main reason for the excessive drainflow recorded in plot 39. Amounts of inflow were estimated by comparing recorded drainflow in plot 39 with the average measured for the remaining 21 plots at the Mellby field. Total inflow in the winter of 1994/5

M.H. Larsson et al. / Journal of Hydrology 215 (1999) 116–134

Fig. 3. A comparison of drainflow from plot 39 with 21 neighbouring plots at Mellby.

was estimated at ca. 165 mm, with periods of zero inflow alternating with influxes of up to 0.15 mm h ⫺1 during peak discharges (Fig. 3) This information was used in the subsequent modelling to specify the bottom boundary condition for water flow. Three groundwater tubes, perforated along their entire length, were also installed, one between plots 39 and 40 and two others outside the plot boundaries. The depth to the water table was recorded manually weekly. Water samples were taken from the groundwater tubes manually weekly, or every second week. All water samples were sent to the laboratory within two weeks, where they were stored at ⫹ 4⬚C until analysed. 2.5. Dye tracing experiment Tracer studies are, to some extent, inconclusive because of their inability to detect spatial patterns of solute movement. A dye tracing study was therefore performed in plot 40 to provide visible evidence for the presence of preferential flow paths at Mellby and to obtain some information on the nature and scale of this process. Ghodrati and Jury (1990) used Acid-Red 1, a red coloured anionic dye, to characterize preferential flow of water in a loamy sand/coarse sand soil which lacked any visible soil structure. They chose Acid-Red 1 after examining the transport behaviour

121

of a number of dyes in the laboratory, because of its high visibility, minimal adsorption and high water solubility. We also used this dye to investigate preferential flow paths in the non-structured sandy soil at Mellby. An area of 6 by 1 m was marked within plot 40 in May 1992. Rain gauges were placed on each corner of the marked plot to measure application rates. Four replicate soil samples (at 5 cm depth intervals to 25 cm depth) were taken to determine the initial water contents. These were very similar to the water content in the surface soil layers at the time of Br ⫺ application in autumn 1994 (ca. 0.21–0.24 m 3 m ⫺3). A 2.67 mm pulse of Acid-Red 1 solution (12.5 g l ⫺1) was then sprayed uniformly over the 6 m 2 area with the help of a pesticide sprinkler. A second sprinkler irrigation using clear tap water, was applied about 14 h after the dye application, and consisted of two separate applications with a time interval between these irrigations of half an hour. The first application lasted for 6 h at an intensity of approximately 2.7 mm h ⫺1. The second application lasted for 3.5 h at a rate of about 6.3 mm h ⫺1. Thus, around 40 mm of irrigation water was added in total to the plot. The irrigation intensity was kept as low as possible in order to simulate normal precipitation rates. Three days after application, a trench was dug to investigate the horizontal and vertical distribution of the dye. Initially, the red colour of the dye did not contrast very well with the brownish black colour of the top soil. Soil evaporation concentrated the dye at the surface, so the dye pattern was more clearly visible after waiting half an hour. Nevertheless, a higher application dose is certainly recommended in future. Ghodrati and Jury (1990) used 320 g Acid Red 1 per m 2, which is ten times higher than our dose. The dye patterns on horizontal cross-sections (1 by 0.4 m) at four different depths (10, 14, 16 and 26 cm) were drawn onto plastic sheets, which were subsequently photographed with a digital camera.

3. Modelling The functional description of two-region preferential flow and transport in the soil matrix outlined in the following text has been programmed into the MACRO model (Jarvis and Larsson, 1998). MACRO also

122

M.H. Larsson et al. / Journal of Hydrology 215 (1999) 116–134

allows the user to predict rapid non-equilibrium fluxes of water and solute in a third, interacting, macropore flow domain. However, in the simulations presented here, the macropore flow domain is almost never activated in the unsaturated zone because the hydraulic conductivity of the sandy matrix at Mellby is sufficiently large to prevent the generation of macropore flow. In addition, the diffusion pathlength regulating mass exchange between macropores and matrix was set to 1 mm (nominal zero) in order to ensure complete equilibrium in solute concentrations in the saturated zone where the macropores are filled with water. The MACRO model is a comprehensive model of water flow and solute transport/transformation processes in field soils. Only those aspects of the model which are relevant to this study are described in the following sections. A full description of the model is given in Jarvis and Larsson (1998). 3.1. Soil water flow The total soil water content u t is given as the sum of the ‘mobile’ and ‘immobile’ water contents, u m and u im respectively.

ut ˆ um ⫹ um :

…1†

If we assume equal bulk densities and porosities in the mobile and immobile regions and also that the hydraulic functions derived in the laboratory on fully saturated soil samples reflect the mobile region in the field (van Dam et al., 1996), then:

ut ˆ F ulab ⫹ Sim 1

…2†

where F is the fraction of the soil matrix volume which contains ‘mobile’ water, u lab is the water content estimated using the laboratory water retention function, Sim is the degree of saturation in the immobile region and 1 is the soil matrix porosity. van Dam et al. (1996) suggested that in water-repellent soil, the fractional volume of soil containing mobile water should vary with the soil water pressure head, decreasing as the soil becomes drier. In this study, the soil water status remained at quasi-steady state throughout the experimental period, so that for the sake of simplicity, F is treated as constant. Both van Dam et al. (1990; 1996) and Saxena et al. (1994) effectively assumed that the degree of saturation in the immobile

region Sim was zero. Later, van Dam et al. (1997) extended the approach to allow a constant non-zero degree of saturation in the immobile region. However, one potential difficulty with this is that if Sim is given a sufficiently large value, a conceptually unrealistic situation can arise in which the total water content becomes larger in the presence of preferential flow, instead of smaller. We therefore introduce the additional assumption that the degree of saturation in the immobile region is a constant fraction k of that in the mobile region, such that: Sim ˆ …1 ⫺ F †k

ulab 1

…3†

and

  ut ˆ F ulab ⫹ …1 ⫺ F†kulab

…4†

or

    1⫺F ut ˆ um 1 ⫹ k : F

…5†

Thus, as we treat both F and k as constants in this study, it follows that the mobile water content is always given as a fixed fraction of the total water content. Richards equation is used to calculate unsaturated flow of ‘mobile’ water in the micropores, with a sink/ source term added to account for the exchange of water between the mobile and immobile regions:      2um 2 2c 2uim ˆ⫺ K ⫹1 ⫺ …6† 2z 2z 2t 2t or        1⫺F 2um 2 2c ˆ⫺ K ⫹1 k 1⫹ F 2z 2z 2t

…7†

where c is the soil water pressure head, z is depth, t is time and K is the unsaturated hydraulic conductivity. In this study, the soil water retention function u lab(c ) is given by the Brooks and Corey (1964) equation:   u ⫺ ur ⫺…1=l† c ˆ cb lab …8† ub ⫺ ur where l is the pore size distribution index, c b is the pressure head at the boundary between macropores and micropores, u b is the equivalent water content and u r is the residual water content. The micropore unsaturated hydraulic conductivity function is given

M.H. Larsson et al. / Journal of Hydrology 215 (1999) 116–134

by Mualem (1976):   u1ab ⫺ ur n⫹2⫹…2=l† K ˆ FKb ub ⫺ ur

…9†

where Kb is the saturated hydraulic conductivity of the micropores, and n is the tortuosity factor. In the macropore region, hydraulic conductivity is given by a simple power law function (Jarvis and Larsson, 1998). An upward flux boundary condition was used at the base of the soil profile in order to solve Eq. (6). This varied during the simulation period from zero to a maximum value of 0.15 mm h ⫺1 depending on the estimated inflow of water to plot 39 from the surrounding area (Fig. 3). 3.2. Soil temperature Soil temperatures are calculated using the heat conduction equation, with the surface boundary condition for snow-free conditions given as the air temperature. With snow cover, the soil surface temperature is calculated assuming steady-state heat flow through a snowpack of constant density and thermal conductivity. Soil temperature at the bottom boundary is calculated from an analytical solution of the heat conduction equation assuming a sinusoidal variation of air temperature on an annual basis. 3.3. Drainage For saturated soil layers above drain depth, the lateral water flow rate qd to primary and secondary drainage systems is calculated using seepage potential theory, accounting for vertical variations in saturated hydraulic conductivity in layered soil (Leeds-Harrison et al., 1986): qd ˆ A f e

…10†

where Af is a geometry or shape factor and e is the seepage potential given by: " # " #! h2u h2l ⫺ Hhl ⫺ …11† e ˆ K Hhu ⫺ 2 2 where H is the water table height and hu and hl are the heights at the top and the bottom of a specific layer, so that hu ⫺ hl equals the layer thickness. Total drainflow is given as the sum of saturated flows from macropores and micropores. These are calculated separately

123

as a function of the saturated hydraulic conductivities K in each respective domain. For a primary drainage system consisting of parallel field drains, the shape factor Af is given by (Youngs, 1980): Af ˆ

8 L2

…12†

where L is the drain spacing. The shape factor for general groundwater seepage to a square-shaped secondary drainage system enclosing a known area Ad is given by (Youngs, 1980): Af ˆ

13:652 : Ad

…13†

Eq. (12) is used for calculating saturated water flow to the primary ‘within-plot’ drains, and Eq. (13) for the secondary ‘protective’ drains surrounding the plot. Drainflow contributions from saturated layers below drain depth are calculated using the second term of the Houghoudt drainage equation. 3.4. Solute transport Solute transport in the ‘mobile’ water is predicted using the convection-dispersion equation written for a non-reactive tracer: ÿ
  X 2 um c m 2 2cm ˆ Dum ⫺ qcm ⫺ U …14† 2z 2t 2z where cm is the solute concentration in the P mobile water, D is the dispersion coefficient and U represents sink/source terms to account for solute losses to primary and secondary drainage systems and lateral groundwater flow and also solute exchange between mobile and immobile regions. The dispersion coefficient represents the effects of both mechanical dispersion and molecular diffusion: D ˆ Dv v ⫹ D0 f*

…15†

where Dv is the dispersivity, v is the pore water velocity ( ˆ q/u m), D0 is the diffusion coefficient in free water and f* is the tortuosity factor calculated according to the model described by Millington and Quirk (1961). Solute losses to primary and secondary field drainage systems are calculated by multiplying the solute concentration in each saturated soil layer with the relevant drainage rate. Thus, it is assumed that the

124

M.H. Larsson et al. / Journal of Hydrology 215 (1999) 116–134

Table 3 Soil hydraulic parameter values derived from measurements Depth (cm)

0–20 20–32 32–60 60–110 a b

Parameter KS a (mm h ⫺1)

Kb (mm m ⫺1)

c b (cm)

u S b (m 3 m ⫺3)

u b (m 3 m ⫺3)

u r (m 3 m ⫺3)

l

66 66 100 100

5 5 14 14

15 15 35 35

0.45 0.42 0.38 0.34

0.39 0.35 0.32 0.32

0.00 0.00 0.05 0.05

0.2 0.2 1.0 1.0

Saturated hydraulic conductivity. Saturated water content.

solute concentration is uniformly mixed laterally within each layer. In addition to leaching losses to field drainage systems, solute loss in lateral shallow groundwater flow Ug is calculated for each saturated soil layer using a retention time concept:   c ⫺ cg …16† Ug ˆ u Rt

mobile and immobile water. The actual field degradation rate is adjusted for soil moisture and soil temperature using the simple functions described by Boesten and van der Linden (1991). A Freundlich sorption isotherm is assumed, with the sorption sites partitioned between mobile and immobile domains according to the relative soil volumes.

where cg is the solute concentration in the local groundwater entering the soil region of interest and Rt is the retention time. The rate of solute exchange US between mobile and immobile regions is given as the sum of diffusive and convective terms:   ÿ
2uim 0 …17† c US ˆ a cim ⫺ cm ⫹ 2t

3.6. Driving variables

where cm and cim are the solute concentrations in the mobile and immobile water, respectively, a is a mass transfer coefficient for solute diffusion and c 0 represents either the macropore or micropore solute concentration depending on the direction of water exchange. In order to solve Eq. (14), boundary conditions at the soil surface and at the base of the soil profile are required. At the soil surface, we assume that applied solute is partitioned between mobile and immobile regions according to the ratio of mobile and immobile water contents at the time of application. A constant solute concentration in the groundwater (cg ˆ 0) is assumed. 3.5. Reactive solute transport For reactive solutes, degradation is calculated using first-order kinetics, with the same half-life assumed in

Measured ‘on-site’ hourly precipitation data was used as input to the model, accounting for wind drift losses using correction factors of 1.07 for rain and 1.14 for snow. Precipitation was assumed to fall as snow when mean daily air temperatures were less than ⫺ 2⬚C and as rain at air temperatures greater than 2⬚C. At intermediate temperatures, the proportion of rain and snow was assumed to vary linearly. Snowmelt was calculated using a simple degree–day concept. Potential evaporation from the bare soil surface during the experiment was calculated on a daily basis from the measured meteorological variables using the Penman–Monteith combination equation, assuming a surface resistance of zero. Actual soil evaporation is estimated from a comparison of this potential rate with the maximum possible water flow rate towards the soil surface (Jarvis and Larsson, 1998). 3.7. Modelling strategy and parameterization The soil hydraulic parameters (see Table 3) were derived by simultaneously fitting the Brooks-Corey and Mualem functions (Brooks and Corey, 1964; Mualem, 1976) to the measured data (Figs. 1 and 2) using the program RETC (Yates et al., 1992). The

M.H. Larsson et al. / Journal of Hydrology 215 (1999) 116–134

125

behaviour in the field, the model was then run with all water as ‘mobile’ (by setting F ˆ 1.0) in order to quantify the impact of preferential flow and transport on solute leaching. A hypothetical pesticide was also simulated, with a degradation half-life at 15⬚C of 10 days in the topsoil (0–32 cm) and 100 days in the subsoil, and a linear isotherm with sorption constants of 0.5 and 0.1 cm 3 g ⫺1 in topsoil and subsoil respectively. The compound was assumed to be applied on the same day as the Br ⫺, while the remaining parameters describing the moisture and temperature response functions were set to default values in the MACRO model. 3.8. Model evaluation In addition to graphical displays of simulations and measured data, an objective statistical measure, the model efficiency (ME), is used to compare simulations with and without preferential flow, and to evaluate the performance of the model, as suggested by Loague and Green (1991): N P

ME ˆ

 2⫺ …O ⫺ O†

1

N P

…P ⫺ O†2

1 N P

 2 …O ⫺ O†

…18†

1

Fig. 4. Observed dye tracing patterns. The dark areas indicate a strong dye colour, while speckled areas indicate non-dyed regions presumed to represent central or ‘core’ regions of fingers.

objective of the model application was to quantify the impact of preferential flow on solute transport. Therefore, we first calibrated the mobile-immobile region model against the field measurements, varying those parameters in the functional model which are difficult to measure or for which no measurements were available (e.g. the mass transfer coefficient a , the fractional soil volume containing mobile water F, the ratio between the degrees of saturation in the immobile and mobile regions k, and the solute retention time in the saturated zone, Rt). The remaining model parameter values were estimated from a combination of literature values (e.g. the diffusion coefficient) or default values in the MACRO model (e.g. dispersivity ˆ 1 cm). Having achieved a satisfactory fit of the mobile–immobile region model to the observed

where P and O represent predicted and observed values respectively and O is the mean of the observed values. The maximum and ideal value for ME is 1.0, while a negative value of ME indicates a poor fit to the data, such that model predictions are worse than using the observed mean as an estimate of the measured data points. 4. Results and discussion 4.1. Evidence of preferential flow 4.1.1. Dye tracing experiment Despite the low irrigation intensity used in the dye experiment (approximately 2.7 mm h ⫺1), the irrigation water immediately started ponding at the soil surface, forming pools ca. 10–20 cm in diameter in small topographical depressions at the soil surface. This is remarkable because the nominal saturated hydraulic conductivity of the topsoil (Fig. 2a) is

126

M.H. Larsson et al. / Journal of Hydrology 215 (1999) 116–134

Fig. 5. Cumulative frequency distribution of bromide, applied to the soil surface and in soil cores sampled to 90 cm on the first two sampling occasions.

more than 20 times larger than the applied irrigation rate. These small pools of water suddenly disappeared about an hour after the irrigation started, indicating that the initial water repellency at the soil surface had been overcome by the pressure head of water at the soil surface, and/or that changes in the wetting contact angle occurred with time. The observed dye patterns are shown in Fig. 4, which clearly indicates the occurrence of preferential flow paths. The percentage of conducting soil volume was largest at the horizontal cross-section at 10 cm depth (ˆ 38%), but a few stained patches could even be observed at 40 cm in the subsoil (not shown). The observed pattern consists of cylinders or fingers. Most water presumably flowed through the non-stained areas of the soil surrounded by stained rings. At the beginning of the irrigation with clean water, the water will mix with the dye at the soil surface and most of this water/dye mixture will move downwards through the saturated inner core of the fingers. Lateral movement of water and dye to the less saturated outer core of the finger will occur because of a concentration gradient (i.e. diffusion) and a pressure potential gradient (i.e. convection). After some time, all the dye on the soil surface will have penetrated into the soil and clean water will infiltrate the soil, also mainly through the inner core of the fingers, and push the water/dye

front further downwards and outwards. The inner core of the fingers in the uppermost part of the soil will then be free of dye. This pattern can be observed in the uppermost 16 cm of the soil. Below this depth, fingers gradually cease to exist. The inner core at the base of each separate finger is stained (see Fig. 4). This is in agreement with the explanation given, as the dye has only just reached this depth (i.e. less time for lateral movement) while no clean water has infiltrated this deep into the soil. At Mellby, finger flow was apparently generated by water repellency and microtopography directly at the air-dry soil surface, which caused the water to redistribute and collect in distinct surface ‘pools’. Using the terminology and model concept of Ritsema and Dekker (1995), the distribution layer may then be considered to be infinitely thin at Mellby. Our experiment also showed the importance of lateral redistribution of water (and solutes) from the fingers to the surrounding dry soil, caused by pressure potential gradients. The rate of this lateral redistribution is not known exactly, but it must have occurred within three days, the time allowed between application and excavation. 4.1.2. Variability of bromide displacement Fig. 5 shows the spatial variation in the Br ⫺ application, together with that observed for the total amount of bromide in soil cores to 90 cm depth on the first two sampling occasions. The centre of mass of Br ⫺was, on average, at 8 cm depth on the first sampling occasion, and at 32 cm on the second. Fig. 5 shows that the degree of spatial variability in Br ⫺ amounts significantly increased from the day of application (CV ˆ 19%) to the first soil sampling 28 days later (CV ˆ 50%) by which time 28 mm of rain had fallen. Thereafter, the variability increased only slightly, with a CV of 59% found 63 days after application. We believe that this can be attributed to lateral redistribution of the Br ⫺ caused by preferential flow generated at or near the soil surface soon after application. This pattern of Br ⫺ redistribution is also consistent with the observations from the dye tracing experiment. This increase in variability cannot be a result of greater variability in Br ⫺ extraction in soil samples compared to the aluminium trays, because such variability would be random and not consistent

M.H. Larsson et al. / Journal of Hydrology 215 (1999) 116–134

127

4.2. Model application

Table 4 Parameters values used in the mobile–immobile model Depth (cm)

F

a (days ⫺1)

k

u m/ u t

0–32 32–70 70–80 80 ⫹

0.7 0.8 0.9 1.0

0.0 0.0 0.0 0.0

0.7 0.7 0.7 1.0

0.77 0.85 0.93 1.00

with depth in each soil cylinder. Ritsema and Dekker (1995) reported significant lateral redistribution of Br ⫺ at the soil surface with only 1 mm of rain following application.

4.2.1. Calibration Table 4 shows the best-fit parameter values for the two-domain mobile–immobile model. Relatively large values of the fractional mobile water content were obtained (e.g. 77% in the topsoil and 85% in the subsoil to 70 cm depth) derived from a fixed value of k of 0.7 throughout the profile, and F values of 0.7 in the topsoil, 0.8 in the subsoil to 70 cm, increasing to 0.9 in the capillary fringe at 70–80 cm and to unity below 80 cm depth (see Table 4). Across the range of values tested (from zero to 0.1 d ⫺1),

Fig. 6. a–d: Measured (mean plus 95% confidence intervals) and predicted soil water contents.

128

M.H. Larsson et al. / Journal of Hydrology 215 (1999) 116–134

Table 5 Model efficiency for water content in soil simulated in one and two domins Sampling occasion and date

2 1994-10-09 3 1994-11-08 4 1994-12-13 5 1995-03-09 Overall

Model efficiency Two-domain case

One-domain case

0.72 0.88 0.45 0.77 0.73

0.32 0.68 0.42 0.93 0.61

model predictions were very insensitive to the value of the diffusive mass transfer coefficient a . This suggests that mass exchange between the mobile and immobile domains was dominated by convection, and supports the field observations discussed in the preceding section. Accordingly, a was set to zero (Table 4). A relatively short retention time for solute in the saturated zone (5 days) was required in order to match the observed Br ⫺ concentrations and leaching rates to the tile drains. This implies a very strong lateral flow across the plot in the shallow groundwater, probably caused by the failure of the drainage systems in nearby plots discussed earlier. 4.2.2. Soil hydrology Figs. 6a–d show comparisons of measured and predicted soil water contents on four sampling occasions, one immediately prior to solute application

Fig. 7. Measured and predicted water table height.

and three in the 150 day winter period following application (Table 2). Table 5 presents the calculations of model efficiency for predictions of soil water contents on each sampling occasion and for the experiment as a whole. For the most part, the efficiencies are high, indicating that the model captures the shapes of the water content profiles, and that the soil water retention curves measured using the evaporation method in the laboratory at only two depths are a reasonable representation of the water retention properties in the field. On the first two sampling occasions (9th October and 8th November), 2 days before and 28 days after Br ⫺ application, the two-domain model better matches the observed soil water contents (Figs. 6a,b; Table 5). Except in the capillary fringe, the predictions by the one-domain model overestimate soil water contents by up to 0.05 m 3 m ⫺3, and generally lie outside the 95% confidence intervals for the mean of the observations. In contrast, on the last sampling occasion at the end of winter (9th March), the one-domain model accurately predicts soil water contents (ME ˆ 0.93), and is clearly better than the two-domain case (Fig. 6d; Table 5). The third sampling occasion in December represents an intermediate situation, where the results indicate incomplete wetting and preferential flow in the upper subsoil at 40–60 cm depth, whereas the measurements in both the topsoil and deeper subsoil match the one-domain simulation better (Fig. 6c). These results suggest that the degree of preferential flow in the soil profile was strongest in early autumn at the time of Br ⫺ application and decreased slowly but steadily throughout the winter, disappearing by the time of the final sampling on 9th March. Nevertheless, the simplifying assumption of constant F and k values in the model seems to represent reasonably well the ‘average’ state during the experiment, probably because soil water conditions did not vary much with time. A functional dependence of F (and possibly k) on soil water content (or pressure head) would certainly have been required if the simulations had extended into the summer period (van Dam et al., 1996). Figs. 7 and 8 show comparisons of measured water table heights and drainflows respectively, with predictions made by the one- and two-domain model simulations. Fig. 7 shows that the model gave good predictions of the water table position, at least until

M.H. Larsson et al. / Journal of Hydrology 215 (1999) 116–134

129

hydrographs during the autumn. This may be due to overestimates of water content in the subsoil (Figs. 6a,b) and suggests that the extent of preferential flow may be somewhat underestimated by the model during this period. Finally, it can be noted that the differences in drainflow predicted by one- and twodomain simulations are marginal, except for the first flow event in autumn, which occurred about one week prior to bromide application.

Fig. 8. Measured and predicted drainflows.

ca. 110 days after Br ⫺ application (late January/early February), and that preferential flow apparently had little measurable impact on the water table. This is most likely because of the large drainable porosity and saturated hydraulic conductivity in this sand soil at drain depth, so that the water table is relatively stable, responding only slowly to water inputs from the unsaturated zone. The large apparent underestimation of water table height in late winter is probably because of the location of the groundwater tube, that is, mid-way between plots 39 and 40. It seems likely that the measurement of water heights during this period was influenced by the failure of the drainage system in plot 40. Fig. 8 shows that the model also gave good predictions of the drainage response, both with regard to the timing of individual flow peaks and their magnitude. However, the model failed to correctly predict the dynamics of drainflow in late winter, from ca. 110 days after bromide application, corresponding to the period when water table positions were also poorly estimated. This reflects the uncertainty concerning the extent of water inflow to the plot. There is also a slight tendency for the model to underpredict drain flow in the first three major

4.2.3. Bromide leaching to drains and mass balance Figs. 9a–c show comparisons of measured and predicted Br ⫺ contents in the soil profile on three sampling occasions following application, while Table 6 presents the calculations of model efficiency for each sampling occasion and for the entire period. The model was able to accurately match the observed mean profile in the surface 20 cm measured on 8th November, 28 days after application (ME ⬎ 0.9, Table 6), although the extent of downward movement of Br ⫺ below 20 cm depth was apparently overpredicted. (Fig. 9a). On this date, the differences between one- and two- domain simulations were small and insignificant in relation to the variability in the measurements (Fig. 9a, Table 6). Differences in the one- and two-domain simulations were more apparent on 13th December, 63 days after Br ⫺ application, when the centre of mass of the observed Br ⫺ pulse had reached ca. 30–40 cm (Fig. 9b). Here, the two-domain model better matched the position of the centre of mass, while the one-domain model underpredicted movement by ca. 10 cm, resulting in a poor (negative) model efficiency (Table 6). Both models overpredicted movement of the leading edge deep in the profile, and over-emphasized tailing of the pulse near the soil surface. It is possible that a smaller value of the dispersivity than the default value of 1 cm would have resulted in improved model predictions. This is also apparent in the profile measured on 9th March, 149 days after application, when the model overestimates Br ⫺ contents in the soil to 60 cm depth, but appears to underestimate contents deeper in the soil profile (Fig. 9c). Fig. 10 shows a comparison of measured and predicted concentrations of Br⫺ in the drainage water. The two-domain model matches the timing of the breakthrough of Br ⫺ to the drains, while the one-domain model predicts a significantly delayed response

130

M.H. Larsson et al. / Journal of Hydrology 215 (1999) 116–134

Fig. 9. a–c: Measured (mean, 95% confidence intervals) and predicted bromide contents in soil.

Table 6 Model efficiency for bromide in soil simulated in one and two domains Sampling occasion and date

3 1994-11-08 4 1994-12-13 5 1995-03-09 Overall

Model efficiency Two-domain case

One-domain case

0.93 0.30 0.31 0.85

0.96 ⫺0.34 0.26 0.81

compared to the measurements. Both models overpredict Br ⫺ concentrations in late-January, ca. 90–110 days after application. In particular, the measurements indicate a marked dilution effect around day 105 (January 23rd 1995), which is not matched by the model, corresponding to a flow peak of ca. 0.4 mm h ⫺1. Measurements of soil hydrology, air and soil temperatures (not presented here), and Br ⫺ concentrations in the drainflow in January 1995 suggest that this dilution effect is caused by preferential flow occurring during snowmelt following a prolonged period of soil freezing (see Johnsson and

M.H. Larsson et al. / Journal of Hydrology 215 (1999) 116–134

Fig. 10. Measured and predicted bromide concentrations in the drainflow.

Lundin, 1991). On freezing, Br ⫺ will be redistributed from the ice-filled larger pores to the smaller pores which still contain liquid water. Precipitation and snowmelt infiltrating the partially frozen soil will move through the large, initially air-filled, pores without coming into contact with the Br ⫺ residing in the smaller pores, thereby causing dilution in the drainflow. This kind of preferential flow cannot yet be simulated with MACRO, as water and heat flow are

131

not coupled in the model. A physically-based twodomain approach for preferential flow in partiallyfrozen soil was described by Sta¨hli et al. (1996). Fig. 11 compares measured and predicted leaching of Br ⫺ to the primary tile drainage system in plot 39, and demonstrates that the two-domain model closely matches the measured loss. In the one-domain case, leaching is somewhat underestimated, mainly because of the failure to match the initial breakthrough to the drains 60–100 days after application. Nevertheless, the effects of preferential flow on Br ⫺ leaching appear marginal in this experiment, reflecting the relatively large fractional transport volume (see Table 4). This is also illustrated in Table 7 which presents the simulated Br ⫺ mass balance, expressed as percentages of the applied amount. Table 7 shows that ca. 46% of the dose was lost to the primary drains by the end of the experimental period on 31st March 1995, and that the two-domain model predicted only ca. 5% more leaching than the one-domain case. The model also estimated that ca. 16% of the applied amount of Br ⫺ was lost in shallow lateral groundwater flow, while leaching to the secondary drains and to deep groundwater were only minor components of the mass balance. This is confirmed by the measurements in groundwater tubes at 1.55 m depth, which showed only marginal increases of Br ⫺ concentrations during the experiment, from 0.97 g m ⫺3 prior to application to a maximum value of 1.09 g m ⫺3 in late February 1995. 4.2.4. Scenario simulations for a reactive solute Table 8 shows the predicted mass balances for the hypothetical pesticide applied at the same time as the Br ⫺. Here, the effects of preferential flow are clearly more pronounced, with the two-domain simulation giving ca. 80% more total leaching than the onedomain simulation. Again, the loss in shallow groundwater flow is a significant proportion of the total loss, which amounts to 1.2% and 2.2% of the application in the one- and two-domain cases respectively.

5. Conclusions

Fig. 11. Measured and predicted bromide loss in tile drainage.

A dye tracing experiment demonstrated that preferential flow occurred in this sandy soil, generated at the soil surface owing to water repellency and microtopography. Additional direct evidence was

132

M.H. Larsson et al. / Journal of Hydrology 215 (1999) 116–134

Table 7 Bromide mass balance, expressed as percentages of the applied amount, at the end of the experiment on 31st March 1995 Component

Measured

Predicted (two-domain)

Predicted (one-domain)

Lost to primary drains Lost to secondary (protective) drains Lost to lateral groundwater flow Leaching to deep groundwater Remaining in soil profile on 9th March a,b Remaining in soil profile a,c

46.4 —

44.7 1.0

42.3 1.0

— — 23.2

16.5 1.2 28.8

15.7 1.1 31.0

—

36.6

40.0

a Corrected for an estimated 0.06 g m ⫺2 bromide in the profile at the start of the experiment, originating from an earlier tracer experiment in 1992. b Measured and predicted to 90 cm depth. c Predicted on 31st March 1995 to 110 cm depth.

Table 8 Mass balance of a hypothetical pesticide, expressed as percentages of the applied amount Component Lost to primary drains Lost to secondary (protective) drains Lost to lateral groundwater flow Leaching to deep groundwater Remaining in soil profile Degraded

One-domain simulation

Two-domain simulation

0.82 0.02

1.53 0.04

0.37 0.004 12.56 86.19

0.65 0.007 12.64 85.10

provided by large increases in the variability of Br ⫺ contents in the soil soon after autumn application. Comparisons between one- and two-domain simulations (i.e., Richards’ equation coupled to the convection-dispersion equation, with and without mobile/ immobile water) clearly showed that preferential flow reduced soil water contents and increased drainflow in the early autumn, leading to a more rapid bromide transport to tile drains. Nevertheless, the quantitative effect of preferential finger flow on total bromide leaching was minor in this experiment, as can also be deduced from the large values of F and u m/u t derived by calibration. Presumably, this is because the solute application was made to relatively wet soil. However, scenario simulations using the calibrated model showed more significant effects for an autumn application of a hypothetical short half-life pesticide, with leaching predicted to increase by ca. 80% as a result of preferential flow. Preferential flow of water

and Br ⫺ was apparently negligible later in the experiment when the soil became wetter. This indicates that the antecedent soil water content may have an influence on the degree of preferential flow and transport, so that the timing of pesticide application will affect the total amount leached, with a higher risk for leaching in initially drier soils.

Acknowledgements We would like to thank the staff at Husha˚llningssa¨llskapet in Halland, especially Magnus Ha˚kansson and Erik Ekre for help in the field, Marianne Cervenka (Department of Ecology and Environmental Research, SLU), Anna Hansson and Amanuel Tesfamikael (Department of Soil Sciences, SLU) for carrying out soil extractions and Br ⫺ analyses, and Jos van Dam (Wageningen Agricultural University) for

M.H. Larsson et al. / Journal of Hydrology 215 (1999) 116–134

helpful discussions during the course of this work. The research reported in this paper was funded by the European Commission – Environment Research Program (Contract EV5V-CT94-0467), within the project ‘Analysis and improvement of existing models of field-scale solute transport through the vadose zone of differently textured soils, with special reference to preferential flow’.

References Bergstro¨m, L., Jarvis, N., Stenstro¨m, J., 1994. Pesticide leaching data to validate simulation models for registration purposes. J. Environ. Sci. Health A29 (6), 1073–1104. Beven, K., Germann, P.F., 1982. Macropores and water flow in soils. Water Resour. Res. 18, 1311–1325. Boesten, J.J.T.I., van der Linden, A.M.A., 1991. Modelling the influence of sorption and transformation on pesticide leaching and persistence. J. Environ. Qual. 20, 425–435. Bond, R.D., 1964. The influence of the microflora on the physical properties of soils. II: Field studies on water repellent sands. Austr. J. of Soil Res. 2, 123–131. Brooks, R.H., Corey, A.T., 1964. Hydraulic properties of porous media. Hydrology paper no. 3, Colorado State University, Fort Collins, Colorado, 27pp. Dekker, L.W., Jungerius, P.D., 1990. Water repellency in dunes with special reference to the Netherlands. Catena Supplement 18, 173–183. Dekker, L.W., Ritsema, C.J., Oostindie, K., Wendroth, O., Pohl, W., Jarvis, N.J., Larsson, M.H., Gaudet, J.-P., 1999. Moisture distributions and wetting rates of soils at experimental fields in the Netherlands, France, Sweden and Germany. J. Hydrol, 215, 4–22. de Rooij, G.H., 1996. Preferential flow in water-repellent sandy soils – Model development and lysimeter experiments. Doctoral thesis, Wageningen Agricultural University, The Netherlands. 229pp. Gamerdinger, A.P., Wagenet, R.J., van Genuchten, M.Th., 1990. Application of two-site/two-region models for studying simultaneous nonequilibrium transport and degradation of pesticides. Soil Sci. Soc. Am. J. 54, 957–963. Ghodrati, M., Jury, W.A., 1990. A field study using dyes to characterize preferential flow of water. Soil Sci. Soc. Am. J. 54, 1558–1563. Glass, R.J., Nicholl, M.J., 1996. Physics of gravity fingering of immiscible fluids within porous media: an overview of current understanding and selected complicating factors. Geoderma 70, 133–164. Glass, R.J., Steenhuis, T.S., Parlange, J.-Y., 1989. Mechanism for finger persistance in homogeneous, unsaturated, porous media: theory and verification. Soil Sci. 148, 60–70. Glass, R.J., Cann, S., Bailey, N., Parlange, J.-Y., Steenhuis, T.S., 1990. Wetting front instability in unsaturated porous media: a

133

three-dimensional study in initially dry sand. Trans. in Porous Media 5, 247–268. Hendrickx, J.M.H., Yao, T.M., 1996. Prediction of wetting front stability in dry field soils using soil and precipitation data. Geoderma 70, 265–280. Hendrickx, J.M.H., Ritsema, C.J., Boersma, O.H., Dekker, L.W., Hamminga, W., van der Kolk, J.W.H., 1991. Motor-driven portable soil core sampler for volumetric sampling. Soil Sci. Soc. Am. J. 59, 27–34. Hendrickx, J.M.H., Dekker, L.W., Boersma, O.H., 1993. Unstable wetting fronts in water-repellent field soils. J. Environ. Qual. 22, 109–118. Hill, D.E., Parlange, J.-Y., 1972. Wetting front instability in layered soils. Soil Sci. Soc. Am. Proc. 36, 697–702. Jarvis, N.J., 1998. Modeling the impact of preferential flow on Nonpoint source pollution. In: Selim, H.M., Ma, L. (Eds.). Physical Nonequilibrium in Soils: Modeling and Application, Ann Arbor Press, Chelsea, Michigan, pp. 195–221. Jarvis, N.J., Messing, I., 1995. Near-saturated hydraulic conductivity in soils of contrasting texture measured by tension infiltrometers. Soil Sci. Soc. Am. J. 59, 27–34. Jarvis, N.J., Larsson, M.H., 1998. The MACRO model (Version 4.1). Technical description http://130.238.110.134:80/bgf/ macrohtm/macro.htm, Department of Soil Science, SLU, Uppsala. Jarvis, N.J., Bergstro¨m, L.F., Brown, C.D., 1995. Pesticide leaching models and their use for management purposes. In: Roberts, T.R., Kearney, P.C. (Eds.). Environmental Behaviour of Agrochemicals, John Wiley and Sons, pp. 185–220. Johnsson, H., Lundin, L.-C., 1991. Surface runoff and soil water percolation as affected by snow and soil frost. J. Hydrol. 122, 141–159. Jury, W.A., Flu¨hler, H., 1992. Transport of chemicals through soil: mechanisms models and field applications. Adv. Agron. 47, 141–201. Kung, K.-J.S., 1990. Preferential flow in a sandy vadose zone: 1. Field observations. Geoderma 46, 59–71. Leeds-Harrison, P.B., Shipway, C.J.P., Jarvis, N.J., Youngs, E.G., 1986. The influence of soil macroporosity on water retention transmission and drainage in a clay soil. Soil Use Manag. 2, 47– 50. Loague, K., Green, R.E., 1991. Statistical and graphical methods for evaluating solute transport models: overview and application. J. Contam. Hydrol. 7, 51–73. McGhie, D.A., Posner, A.M., 1980. Water repellency of a heavytextured western Australian surface soil. Austr. J. Soil Res. 18, 309–323. Millington, R.J., Quirk, J.P., 1961. Permeability of porous solids. Trans. Faraday Soc. 57, 1200–1207. Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12, 513–522. Nieber, J.L., 1996. Modeling finger development and persistence in initially dry porous media. Geoderma 70, 207–229. Richardson, J.L., 1984. Field observations and measurements of water repellency for soil surveyors. Soil Survey Horiz. 25, 32–36.

134

M.H. Larsson et al. / Journal of Hydrology 215 (1999) 116–134

Ritsema, C.J., Dekker, L.W., 1995. Distribution flow: A general process in the top layer of water repellent soils. Water Resour. Res. 31, 1187–1200. Saxena, R.K., Jarvis, N.J., Bergstro¨m, L., 1994. Interpreting nonsteady state tracer breakthrough experiments in sand and clay soils using a dual-porosity model. J. Hydrol. 162, 279–298. Sta¨hli, M., Jansson, P.-E., Lundin, L.-C., 1996. Preferential water flow in a frozen soil – a two-domain model approach. Hydrol. Process. 10, 1305–1316. Stolte, J., Veerman, G.J., Wopereis, M.C.S., 1992. Manual soil physical measurements, version 2.0, Technical Document 2, DLO Winand Staring Centre, Wageningen, Netherlands, 49pp. van Dam, J.C., Hendrickx, J.M.H., van Ommen, H.C., Bannink, M.H., van Genuchten, M.Th., Dekker, L.W., 1990. Water and solute movement in a coarse-textured water-repellent field soil. J. Hydrol. 120, 359–379. van Dam, J.C., Wosten, J.H.M., Nemes, A., 1996. Unsaturated soil water movement in hysteretic and water repellent field soils. J. Hydrol. 184, 153–173. van Dam, J.C., Huygen, J., Wesseling, J.G., Feddes, R.A., Kabbat, P., van Walsum, P.E.V., 1997. SWAP User’s manual. Simula-

tion of transport processes in the Soil-Water-Air-Plant environment. Report, Dept. of Water Resources, Wageningen Agricultural University, Wageningen, Netherlands. van Genuchten, M.Th., Wierenga, P., 1976. Mass transfer studies in sorbing porous media. I: Analytical solutions. Soil. Sci. Soc. Am. J. 40, 473–480. van Genuchten, M.Th., Wagenet, R.J., 1989. Two-site/two-region models for pesticide transport and degradation: Theoretical development and analytical solutions. Soil Sci. Soc. Am. J. 53, 1303–1310. Wind, G.P., 1966. Capillary conductivity data estimated by a simple method. Water in the unsaturated zone. Symp. 1966, Proc. UNESCO/IASH, 181–191, Wageningen, Netherlands. Yates, S.R., van Genuchten, M.Th., Warrick, A.W., Leij, F.J., 1992. Analysis of measured, predicted and estimated hydraulic conductivity using the RETC computer program. Soil Sci. Soc. Am. J. 56, 343–354. Youngs, E.G., 1980. The analysis of groundwater seepage in heterogeneous aquifers. Hydrol. Sci. Bull. 25, 155–165. Zu¨rmuhl, T., Durner, W., 1996. Modeling transient water and solute transport in a biporous soil. Water Resour. Res. 32, 819–829.