In Denmark more than 95% of the drinking water production

TECHNICAL REPORTS: VADOSE ZONE PROCESSES AND CHEMICAL TRANSPORT Hydrogeological Relationships of Sandy Deposits: Modeling of Two-Dimensional Unsatura...
Author: Mabel Shepherd
3 downloads 1 Views 1MB Size
TECHNICAL REPORTS: VADOSE ZONE PROCESSES AND CHEMICAL TRANSPORT

Hydrogeological Relationships of Sandy Deposits: Modeling of Two-Dimensional Unsaturated Water and Pesticide Transport Bo V. Iversen* University of Aarhus Peter van der Keur and Henrik Vosgerau Geological Survey of Denmark and Greenland Prediction of the movement of water and solutes in the vadose zone requires information on the distribution of spatial trends and heterogeneities in porous media. The present study describes different lithofacies origination mainly from glaciofluvial deposits. Among different lithofacies, hydrological relationships were investigated. By means of a two-dimensional hydrological model it was evaluated how the flow of water and leaching of metribuzin (4-amino-6-tert-butyl-4,5-dihydro-3-methylthio1,2,4-triazin-5-one) was affected. Two selected large outcrop sections consisting of glacial outwash deposits were used in the modeling study. Eleven different lithofacies were distinguished and described in terms of texture distribution, sorting, bedding style, and external boundaries based on excavated soil profiles from 27 locations representing seven predominantly sandy landforms in Denmark. Undisturbed soil columns were sampled from each of the lithofacies and brought to the laboratory to be analyzed. With respect to their soil hydraulic properties, the different lithofacies formed four different hydrofacies having relatively homogeneous, hydrogeological properties. Two large outcrop sections from one of the locations (a gravel pit) located near the terminal moraine of the former Weichsel glacier were used for the HYDRUS-2D modeling. Modeling results revealed that the spatial distribution of sedimentary bodies affected water flow and the leaching of metribuzin.

Copyright © 2008 by the American Society of Agronomy, Crop Science Society of America, and Soil Science Society of America. 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. Published in J. Environ. Qual. 37:1909–1917 (2008). doi:10.2134/jeq2006.0200 Received 22 May 2006. *Corresponding author ([email protected]). © ASA, CSSA, SSSA 677 S. Segoe Rd., Madison, WI 53711 USA

I

n Denmark more than 95% of the drinking water production is based on groundwater. Pesticides and their degradation products during the past decade have increasingly been detected and are now present in much of the groundwater (Kjær et al., 2005a). In 2000, the Danish Ministry of the Environment initiated a large project with the overall objective of providing the necessary background knowledge for identifying areas particularly prone to pesticide leaching (Nygaard et al., 2004). The focus is on sandy agricultural areas that cover about 58% of the land area in Denmark, where approximately 11% of the area in Denmark consists of glaciofluvial deposits (Madsen et al., 1992). Glacial outwash rivers are typically of multiple-channel (braided) type and are normally of low sinuosity (Miall, 1983). Sediments from glaciofluvial deposits are composed of stratified gravel and sand and are a reflection of the dominance of braided stream activity in this environment. Grain sizes decrease from gravel-dominated near the ice front to sand or silt in distal reaches. Hydrogeologists often focus on the quantification of the heterogeneity of the material where the water is flowing. Much of the sedimentological literature, however, consists of qualitative descriptions of sedimentary units and qualitative approaches to describing units, while hydrogeologists require quantitative descriptions (Anderson et al., 1999). Hydrogeologists focus to a high degree on the quantification of different geological materials with respect to their hydraulic properties. When describing glaciofluvial deposits, the sediment can be divided into different lithofacies according to the texture distribution, sorting, bedding style (arrangement of geological layers), and external boundaries. Relating lithofacies to textural and structural units with unique hydraulic properties will provide the hydrogeologist with information on the distribution and heterogeneities of the porous media. The concept of relating lithofacies in gravel deposits to different hydrogeological facies was described by Anderson (1989). A hydrogeologic facies (hydrofacies) is defined to be a homogeneous, but anisotropic unit that is hydrogeologically meaningful for purposes of field experiments and modeling (Anderson, 1989). The use of hydrofacies, therefore, allows hydrogeologists to predict the expected B.V. Iversen, Univ. of Aarhus, Faculty of Agricultural Sciences, Dep. of Agroecology and Environment, Research Centre Foulum, PO Box 50, DK-8830 Tjele, Denmark. P. van der Keur and H. Vosgerau, Geological Survey of Denmark and Greenland, Dep. of Hydrology, Øster Voldgade 10, DK-1350 Copenhagen, Denmark. H. Vosgerau, current address, Miljøcenter Roskilde, Ny Østergade 7-11, DK-4000 Roskilde, Denmark. Abbreviations: k, unsaturated hydraulic conductivity; Ks, saturated hydraulic conductivity; CV, coefficient of variation.

1909

Table 1. Description of the studied sites. Site Mjøls Ragners gravel pit

Deposit glaciofluvial glaciofluvial

Ap and B† x

Described facies¶ Gd, St Gd, Gs, Gt, Gv, Sl, St, Sv‡

Stubkær glaciofluvial x Sp Knivsig glaciofluvial St Kølvrå glaciofluvial Sp Røjen glaciofluvial Gd, St, Sv Emmerske glaciofluvial Sh Hallundbæk§ glaciofluvial St Hogager glaciofluvial St Ommose§ glaciofluvial Sr Simmelkær glaciofluvial Sp Frederiks glaciofluvial x – Nedre Julianehede§ glaciofluvial x – Nørlund§ glaciofluvial x – Søbjerg§ glaciofluvial x – Ruskær§ glaciofluvial x – Alle glaciofluvial Sh Hjortkær-DS glaciofluvial Sv Nr. Felding-DS glaciofluvial Sh Hjortkær-MS till Dmh Nr. Felding-MS till Dmh Gammelkirke till Dmh Låstrup till Dmh Sjørup till Dmh Astrup§ till Sh Ajstrup marine Sh Ulsted marine Sr, Sv † Hydraulic properties measured in the Ap and B horizon. ‡ Hydraulic properties not measured in Gs, Gt, Gv, and Sl. § Hydraulic properties measured on 100 cm3 soil cores only ¶ Gd, Gravel (homogeneous or poorly layered); Gs, Gravel (planar x-beds); Gt, Gravel (trough x-beds); Gv, Gravel (horizontal bedded); Sr, Sand (ripple x-laminated); Dmh, Diamicts (matrix supported, homogeneous); Sh, Sand (homogenous); Sl, Sand (x-beds, low angle); Sp, Sand (planar x-beds); St, Sand (trough x-beds); Sv, Sand (horizontal laminated).

spatial trends in hydraulic conductivity from limited site-specific information. With the obtained knowledge of large-scale trends in sediment facies, these trends can be matched with site-specific information on hydraulic conductivity if it is possible to isolate a sediment body that is homogeneous with respect to hydraulic conductivity. Until today the study of the relationship of lithofacies with hydrofacies has focused on groundwater flow and groundwater contamination only and the characterization of the hydraulic properties has focused on measuring the saturated hydraulic conductivity only (Anderson et al., 1999; Bersezio et al., 1999; Klingbeil et al., 1999; Heinz et al., 2003). Little attention has been paid to the hydraulic properties measured in the vadose zone. The flow of water in the vadose zone is strongly affected by the heterogeneity of the soil. Several investigations have revealed that preferential flow paths constitute the dominant flow in heterogeneous sandy soils where layers with coarse sand and layers with more fine-textured materials frequently alternate (Kung,

1910

1990; Walter et al., 2000; Alfnes et al., 2004). On sloping layered soil where this kind of textural change is seen (fine-textured over coarse-textured layers), water is forced to flow laterally on top of the underlying coarse layer resulting in a capillary barrier. This type of preferential unsaturated flow in the vadose zone is commonly referred to as funneled flow (Kung, 1990) and might be a major cause of the contamination of groundwater. Ju and Kung (1997) performed numerical simulations with a finite element model (Ju and Kung, 1993). They studied the impact of funnel-type preferential flow on contaminant transport in four 12-m wide by 6-m deep two-dimensional hypothetical profiles with embedded coarse sand layers and found that contaminant transport (bromide) along funnel-type flow paths was significantly accelerated. They also discovered that the impact of funnel flow on contaminant transport was most severe when the net infiltration was low. In another experiment (Alfnes et al., 2004), numerical simulations of a field experiment on a glacial delta with sloping layers of alternating finer and coarser sand were performed. The results showed only a small capillary barrier effect in the simulations at the interface between a fine sand and a coarse sand layer. The study indicated that anisotropy was more important for the heterogeneity of the flow than capillary barriers even at low flow rates. The present study aims, first, to describe different lithofacies from soils originating mainly from glaciofluvial deposits and to investigate hydrogeological relationships among the different lithofacies. The focus will be on the unsaturated hydraulic properties of the deposits. Second, the response characteristic in the flow of water and leaching of pesticides in the vadose zone in two selected large outcrop sections consisting of glacial outwash deposits will be evaluated by means of a two-dimensional hydrological model.

Materials and Methods Soils Studied In the study, 27 different agricultural sandy soils representing seven predominantly sandy landforms (geomorphological units) in Denmark were described. Nineteen soils originated from glaciofluvial deposits, six soils originated from till deposits, and two soils originated from a marine deposit (Table 1). In 22 of the profiles, 11 different lithofacies were distinguished and described on the basis of texture distribution, sorting, bedding style, and external boundaries. The description of the lithofacies was based on Miall (1977). A description of the investigated lithofacies is given in Table 2. Soil samples from each described soil horizon and lithofacies were brought to the laboratory to be analyzed for soil textural distribution and amount of organic matter. For most lithofacies, five large undisturbed 6.28 × 10−3 m3 soil cores (0.20-m diam., 0.20-m length) and five small 10−4 m3 (0.061-m diam., 0.034-m length) soil cores were extracted. For extracting the 6.28 × 10−3 m3 cores, steel columns were gently forced into the soil with a hydraulic press mounted on a tractor. For the 10−4 m3 cores, steel rings were forced into the soil by means of a hammer, using a special flange. A few of the 6.28 × 10−3 m3 steel rings were forced into the soil also by using a hammer and a flange, since it was not possible to use the tractor. The depth of the observed facies from where soil samples where extracted varied between 0.40 and 2 m measured from the soil

Journal of Environmental Quality • Volume 37 • September–October 2008

Table 2. Description of the studied lithofacies. Facies†

Description

Gd Gravel, homogeneous or poorly layered Gs‡ Gravel, planar x-beds§ Gt‡ Gravel, trough x-beds Gv‡ Gravel, horizontal bedded Sr Sand, ripple x-laminated¶ Dmh Diamicts, matrix supported, homogeneous#†† Sh Sand, homogenous Sl‡ Sand, x-beds, low angle Sp Sand, planar x-beds St Sand, trough x-beds Sv Sand, horizontal laminated † Based on Miall (1977). ‡ Hydraulic properties not measured. § x-beds: cross beds, inclined sedimentary structures thicker than 1 cm. ¶ x-laminated: cross laminated, inclined sedimentary structures thinner than 1 cm. # Diamicts: poorly sorted gravel-sand-mud deposit. †† Matrix supported: fine grains support the large grains.

surface. After careful removal of the soil-filled cylinder, end surfaces were trimmed with a knife. The samples were protected from evaporation and physical disruption and brought to the laboratory where they were stored at 2 to 5°C until analysis. In the laboratory, the large, undisturbed 6.28 × 10−3 m3 soil cores were slowly saturated in a drainage box with a ceramic plate for 3 d. After that, the cores were drained to a matric potential of –20 hPa and removed to a drip infiltrometer where they were placed on a ceramic plate. The drip infiltrometer (van den Elsen et al., 1999) measured the unsaturated hydraulic conductivity (k [m s−1]) in the wet (near saturated) range at matric potentials (h [hPa]) between approximately −1 to −100 hPa. Five tensiometers were placed in the soil core and a needle device applied a water flux to the top of the sample. The measurement was started as close to saturation as possible while applying a constant suction to the bottom of the soil core and a constant water flux to the top of the sample. When a steady-state pressure head was reached, the flux density was determined. This, together with measurements of the pressure head gradient, yielded the hydraulic conductivity at the resulting matric potential. Subsequently, the flux was decreased and the procedure described above repeated measuring at a lower pressure head. After the measurement of k, the soil cores were placed in the device used for measurement of saturated hydraulic conductivity (Ks [m s−1]). The soil core was resaturated overnight and Ks was determined using the constant head method suggested by Klute and Dirksen (1986). For further information on the sampling procedure of the soil cores and the measurement of k and Ks consult Iversen et al. (2004). The 10−4 m3 soil cores were weighed and placed on top of a sandbox and saturated with water from beneath. The soil water characteristics were determined by draining the soil sample successively to matric potentials of −10, −16, −50, −100, −160, −500, −1,000, and −16 000 hPa using a sandbox for matric potentials from −10 to −100 hPa and a ceramic plate for matric potentials from −160 to −16,000 hPa. Water content at −16,000 hPa was

Iversen et al.: Hydrogeological Relationships of Sandy Deposits

Fig. 1. Description of Section I from the outcrop section with (a) the distribution of the investigated lithofacies, (b) the resulting hydrofacies, and (c) normalized pesticide concentration at a depth of 3 m (dashed line in b) after 3, 5, and 10 yr of typical metribuzin applications to potatoes intercropped with spring barley. Gd, Gravel (homogeneous or poorly layered); Gt, Gravel (trough xbeds); Sl, Sand (x-beds, low angle); St, Sand (trough x-beds).

determined after the soil had been ground and sieved through a 2mm sieve (Klute, 1986). When Ks had to be measured on the 10−4 m3 soil cores, the soil samples were resaturated and the conductivity was measured using the constant head method (Klute and Dirksen, 1986). Finally, the soil samples were oven-dried at 105°C for 24 h and weighed to determine the bulk density and water content.

Two-Dimensional Modeling The HYDRUS-2D model (ver. 2.05) is a finite element model for simulating the movement of water, heat, and multiple solutes in variably saturated media (Simunek et al., 1999). The program numerically solves the Richards equation for saturated-unsaturated water flow and the advection-dispersion equations for solute transport. The flow equation incorporates a sink term to account for water uptake by plant roots. For the modeling, two large outcrop sections (Section I and II) from one of the locations (Ragners gravel pit) where the sediment was exposed were used (Fig. 1a and 2a). The gravel pit was located near the terminal moraine of the former Weichsel glacier

1911

Fig. 2. Description of Section II from the outcrop section with (a) the distribution of the investigated lithofacies, (b) the resulting hydrofacies, and (c) normalized pesticide concentration at a depth of 3 m (dashed line in b) after 3, 5, and 10 yr of typical metribuzin applications to potatoes intercropped with spring barley. Gd, Gravel (homogeneous or poorly layered); Gs, Gravel (planar x-beds); Gt, Gravel (trough x-beds); Gv, Gravel (horizontal bedded); Sl, Sand (x-beds, low angle); Sr, Sand (ripple x-laminated); St, Sand (trough x-beds).

in Denmark. The outcrop sections consisted of glacial outwash deposits from the former fluvial drainage zone of the glacier. The two outcrop sections were described according to the distribution of recognized lithofacies. To represent the soil water retention and hydraulic conductivity functions, parameters of the van Genuchten (1980) and Mualem (1976) models were fitted to measured data using the RETC code (van Genuchten et al., 1991). The equation of van Genuchten (1980) states that:

m k (θ ) = K s S e1/2 ⎡⎢ 1 − (1− S e1/m ) ⎤⎥ ⎣ ⎦

2

[3]

where θ (m3 m−3) is the water content and θr (m3 m−3) and θs (m3 m−3) are the residual and saturated water contents, respectively. The Mualem-van Genuchten equation describing the hydraulic conductivity is defined as:

The measured hydraulic data were used for the estimation of hydraulic parameters using a relative weight of conductivity data against retention data of 1:5, since it was expected that the conductivity data showed more scatter than the water content data (van Genuchten et al., 1991). The optimization was performed on hydraulic data from each of the investigated facies (see Results and Discussion section). An optimization of hydraulic parameters was also performed on the hydraulic data for the Ap and B horizons from the excavated profiles. As an upper boundary condition, climatic data from the Tylstrup climate station were used with 8.493 m rainfall for the 10-yr simulated period from 1 Jan. 1991 to 31 Dec. 2000. Mean annual rainfall at Tylstrup is similar to the mean annual precipitation averaged over Denmark. To make the simulations as realistic as possible, the hydraulic properties of the Ap and B horizons measured at six sites all located near the terminal moraine of the former Weichsel glacier were used as input to the model. The Ap horizon was given a thickness of 0.30 m and the B horizon a thickness of 0.55 m in the uppermost part of the outcrop section,

1912

Journal of Environmental Quality • Volume 37 • September–October 2008

n S e = ⎡⎢ 1 + (αh ) ⎤⎥ ⎣ ⎦

−m

[1]

where α (m−1), n, and m are empirical parameters. In the present study, m was set to 1 − 1/n. The effective degree of saturation, Se, is defined as: Se =

θ −θ r θ s −θ r

[2]

corresponding to an average of the measured horizons in the field. To initialize the model with respect to the soil water distribution, a “warm-up” period of 1 yr (1990) was used. Agricultural management was defined in the model for a typical metribuzin (4-amino-6-tert-butyl-4,5-dihydro-3-methylthio-1,2,4-triazin5-one) application to potatoes intercropped with spring barley. A total amount of 0.245 kg active metribuzin per hectare was applied to the field in the middle of June each year when growing potatoes. Since HYDRUS-2D does not have a crop-growing module, it cannot subdivide potential evaporation from the climatic data into potential evaporation and potential transpiration. Instead, the model EVACROP (Olesen and Heidmann, 1990) was used to convert the climatic data into the respective variables. Values of the adsorption isotherm coefficient, Kd (Table 3), were measured at the nearby locality at Stubkær (Torp, 2005). The degradation rate was determined by taking in situ soil samples in which degradation is assumed to take place in the dissolved phase where the liquid and adsorbed phase is in equilibrium. Degradation in the adsorbed pool (solid) is assumed negligible and set to zero. Values of the first-order degradation rate constant for the dissolved phase, μw (Table 3), were taken as 80th percentiles because rather low estimated degradation rates tend to be more uncertain and less reliable than the others (van der Keur et al., 2004). A Freundlich adsorption isotherm exponent, β, of 0.9 was used in the modeling. The longitudinal and the transverse dispersivities were set to 0.50 and 5 × 10−4 m, respectively, which corresponds to the findings of Christiansen (2003). Here, the longitudinal dispersivity was experimentally determined by calibrating the advection-dispersion component of HYDRUS-2D against observed bromide content in suction cells for two soil profiles at the Tylstrup location. A value for the longitudinal dispersivity of 0.50 m was experimentally found to best represent the two profiles and is in the same order of magnitude found in other experiments in Denmark under similar conditions and also found by Vereecken et al. (1995) for a German site. This value was subsequently satisfactorily validated by simulation bromide concentration using HYDRUS-2D for another nearby profile at 1-m and 2-m depth and compared with measured values. The estimated dispersivity increases with distance, which is known as the scale dependency, but is not important here for the vadose zone. A transverse dispersivity of 5 × 10−4 m was also found by Jensen et al. (1993). The modeling parameters with respect to the pesticide were used irrespective of the individual hydrofacies since the hydrofacies used for modeling all had low content of clay and organic matter beneath the Ap and B horizons. The groundwater level was not observed 10 m below ground surface at the nearby locality at Frederiks (Torp, 2005). Therefore, lower boundary conditions were set to free drainage where a unit total vertical hydraulic gradient was assumed (gravity flow only). Simulations were not calibrated against measured field data of pesticide leaching, as the main interest was to use the model results in a relative sense to investigate the effect of pesticide leaching when taking into account the heterogeneity of the sandy outcrop sections.

Iversen et al.: Hydrogeological Relationships of Sandy Deposits

Table 3. Modeling parameters related to the properties of metribuzin. Horizon

Kd† 3

μw‡ −1

m kg 0.00343 0.00017 3.33 × 10−5

s−1 5.3 × 10−8 2.1 × 10−9 2.1 × 10−9

Ap B C † Adsorption isotherm coefficient. ‡ First-order degradation rate constant for the dissolved phase.

Results and Discussion Most of the gravel facies (Gd, Gs, Gt, and Gv) were located close to the terminal moraine of the former Weichsel glacier in Denmark on the glacial outwash plain (Table 1). Only one gravel facies was found on the most distal parts of the outwash plain (Røjen). At Stubkær, which also was located near the terminal moraine, no gravel facies were found. For the remaining sites with glaciofluvial deposits, different forms of sandy facies were found (Sr, Sh, Sl, Sp, St, Sv). No diamicts (poorly sorted gravel-sand-mud deposits) were found for the glaciofluvial deposits, but were only related to till deposits. All described facies were found in the C horizon of the described profiles, since the bedding style and structure of any facies in the A or B horizon was obliterated due the soil pedogenetic processes. The soil water characteristics showed (Fig. 3a) that for the gravel facies (Gd) the water content decreased sharply from saturation until a matric potential of −50 hPa where most of the water had drained out of the soil. The sandy facies showed a marked drainage of the soil from a matric potential of −10 hPa until −160 hPa. This marked drainage is a reflection of the sandy/gravelly texture of the soil (Fig. 4) leading to a large number of soil pores having an equivalent pore diameter larger than 20 μm. For the diamicts facies (Dmh) the soil water characteristics showed a more gradual drainage of the soil indicating a more equal distribution of soil pores. For saturated hydraulic conductivity (Ks) the largest deviations were observed for the Dmh facies, whereas the smallest deviations were observed for the Gd facies (Fig. 3b). The Dmh facies had values of Ks ranging from less than 10−7 m s−1 to 5 × 10−4 m s−1. Highest values of Ks were observed for the Gd facies with up to 10−3 m s−1. For the unsaturated hydraulic conductivity (k), the lowest values and largest scatter of conductivities were observed for the Dmh facies, which also had the lowest values measured. Probably, the Dmh facies had the highest number of structural macropores, as suggested by the large difference in the conductivity between the saturated and the near-saturated hydraulic conductivity for some of the soil samples. The large difference in the conductivity was to a large extent explained by flow through macropores when the soil is fully saturated. Low values of Ks in the Dmh facies were probably related to soil samples without larger macropores. The other sandy or gravelly facies showed a steep decline in the hydraulic conductivity between −10 and −100 hPa by almost three orders of magnitude (10−4 to 10−7 m s−1). These facies did not show large differences in the conductivity between the saturated and the near-saturated measurements, which implies that there is limited macropore flow in these soil types when the soil is fully saturated.

1913

For the fitted van Genuchten parameters, low values and large confidence intervals (95%) were obtained for θr (Fig. 5a). For the fitted values of θs, large differences were observed between the Dmh facies and the other facies (Fig. 5b). Smaller differences were observed between the sandy facies. The highest values of θs were observed for the Ap horizon. This is related to the higher content of organic matter leading to a higher soil porosity. The α parameter (defining the inflection point of the soil water characteristic curve) showed high values of facies Gd and low values for facies Dmh, Sh, and Sr (Fig. 5c). For the n parameter (the slope parameter of the soil water characteristic curve) low values were observed for the Dmh facies (Fig. 5d). More or less equal values of n were observed for the sandy and gravelly facies. The plotted van Genuchten functions for the seven different facies are seen in Fig. 6a and 6c. Based on the shape of the soil water retention characteristics and the hydraulic conductivity curve, and the interpretation of the statistics in Fig. 5, the seven lithofacies were divided into four distinct hydrofacies (Fig. 6b and 6d). Especially Hydrofacies 3 (facies Dmh) differed from the three

other hydrofacies with low values of θs, α, and n. Hydrofacies 1 (facies Gd) showed high values of α and high values of k near saturation. Hydrofacies 2 (facies Sr) showed low values of α, whereas Hydrofacies 4 (facies Sh, Sp, St, and Sv) all showed intermediate values of α and k near saturation. As seen from Fig. 4, differences in the hydraulic properties were mainly related to differences in the soil texture distribution and to a lesser extent to differences in the depositional structure of the lithofacies. The high content of stones and coarse sand in facies Gd led to a high content of large soil pores with a low water holding capacity (large quantity of water drained out of the soil from saturation to about −100 hPa). On the other hand, facies Dmh had the highest content of small soil particles and was at the same time the most unsorted facies. This led to a higher content of small soil pores, a more uniform distribution of soil pores, and a higher water holding capacity of the soil. Facies Sr differed from the other sandy facies by having the highest content of fine sand and being the most well sorted. This was reflected in the hydraulic properties of facies Sr with a high water holding capacity compared to the other sandy facies. This relatively high water holding capacity was caused by a high content of intermediary sized pores. Only sandy and gravelly facies were found in the two outcrop sections used for the modeling. In Section I, four different facies (Gd, Gt, Sl, and St) were found (Fig. 1a). In Section II, seven different facies (Gd, Gs, Gt, Gv, Sl, Sr, and St) were found (Fig. 2a). The large number of gravelly facies was in agreement with the location of the outcrop sections near the terminal moraine of the former Weichsel glacier. Here, the facies were typically deposited by braided rivers with a high meltwater discharge. The layers were relatively horizontally distributed in the profile. Some of the gravelly facies had a sloping end section. For the gravelly facies, only the hydraulic properties of Gd were described (Table 1). For the sandy facies, the hydraulic properties of all facies from the gravel pit were described, except for Sl (Table 1). In relation to the modeling, all gravelly facies were given the same hydraulic properties. The sandy facies Sl was given the same hydraulic properties as the other sandy facies related to Hydrofacies 4. Therefore, in Section I (Fig. 1b), two hydrofacies were represented in the profile (Hydrofacies 1 and 4). In section II (Fig. 2b) three different hydrofacies were represented in the profile (Hydrofacies 1, 2, and 4), where Hydrofacies 2 was represented at the bottom of the profile only. Modeled concentrations of metribuzin at a depth of 3 m after a simulation period of 3, 5, and 10 yr are presented in Fig. 1c and 2c. Values are normalized with respect to the highest measured metribuzin concentrations in one of the sections at that depth during the modeling. After a simulation period of 3 yr the highest metribuzin concentrations were found in Section I (Fig. 1c). Here normalized concentrations of metribuzin varied between 0.54 and 0.78, whereas normalized concentrations of metribuzin in Section II (Fig. 2c) varied between 0.42 and 0.70. The average normalized concentrations after 3 yr for Section I and II were 0.67 and 0.52, respectively. For Section I, the lowest normalized con-

1914

Journal of Environmental Quality • Volume 37 • September–October 2008

Fig. 3. (a) Measured soil water characteristic (mean values) and (b) measured saturated (Ks) and unsaturated (k) hydraulic conductivity for investigated lithofacies. Error bars show ± 1 standard deviations. Gd, Gravel (homogeneous or poorly layered); Sr, Sand (ripple x-laminated); Dmh, Diamicts (matrix supported, homogeneous); Sh, Sand (homogenous); Sp, Sand (planar x-beds); St, Sand (trough x-beds); Sv, Sand (horizontal laminated).

centrations after 3 yr were found on the left side of the profile. The sharp increase in the normalized concentration from 0.8 to 1.4 m was probably related to the slightly sloping structure of the two contrasting Hydrofacies 1 and 4, which lead to a preferential flow path due to a small capillary effect. For Section II, peaks in the normalized metribuzin concentrations along the profile were observed at 4.8 and 8.8 m, but it was difficult to relate the normalized concentration peaks to a specific sloping layer. The heterogeneous profile with frequent alternations between Hydrofacies 1 and 4 in the profile was probably responsible for a blurring of the capillary barrier effects in relation to the sloping layers. Overall, only relatively small capillary barrier effects were observed for the two sections. This was related to the relatively small differences in the hydraulic properties beFig. 4. Soil texture distribution of the investigated lithofacies (mean values). Gd, Gravel tween Hydrofacies 1 and 4. Both facies were (homogeneous or poorly layered); Sr, Sand (ripple x-laminated); Dmh, Diamicts (matrix supported, homogeneous); Sh, Sand (homogenous); Sp, Sand (planar x-beds); coarse-textured and the contrast between St, Sand (trough x-beds); Sv, Sand (horizontal laminated). them was too small to observe larger differences in flow patterns. Only at low flow rates 1 μg L−1. The higher simulated metribuzin values in the current when the matric potential was lower than −10 hPa, were there work are probably related to differences in the sediment (coarse differences in the hydraulic conductivity of more than two sand and gravel vs. fine sand) leading to a higher hydraulic conorders of magnitude. At flow rates higher than −10 hPa, only ductivity. In the present study we are focusing on relative values relatively small differences in the hydraulic conductivity were of leached pesticide for different hydrofacies and no calibration observed between the two hydrofacies. Already after 5 yr, difof the model was conducted. Therefore, direct comparison of ferences in the metribuzin normalized concentrations along non-normalized simulated values to threshold values like the 0.1 the 3-m line diminished (Fig. 1c and 2c). After simulation μg L−1 is likely to be misleading. Furthermore, differences are periods of 7 yr, an upper level of metribuzin normalized concentrations was reached. At the end of the simulation period, the mean normalized concentrations between the two sections became the same (Fig. 1c and 2c). This diminishing difference along the profile was probably related to the effect that an increased accumulated amount of leachate at 3-m depth applied over time resulted in an spatial averaging out of the differences in hydraulic properties of the hydrofacies as well as of small-scale heterogeneities at that depth leading to a higher and more evenly distributed normalized concentration profile. The average values of simulated amounts of metribuzin at a depth of 3 m after a simulation period of 10 yr corresponded to non-normalized concentrations of 17 μg L−1. Kjær et al. (2005b) evaluated leaching of metribuzin and its primary metabolites at a fine sandy test site Fig. 5. (a–d) Fitted van Genuchten parameters (θr, θs, α, and n) with 95% confidence in Denmark and found values leached from the intervals (from RETC output) for the investigated lithofacies and the Ap and B root zone in average concentrations considerhorizons. θr (residual water content), θs (saturated water content), α and n are fitting parameters of the Mualem-van Genuchten equation. Dmh, Diamicts (matrix ably exceeding the EU limit value for drinking supported, homogeneous); Gd, Gravel (homogeneous or poorly layered); Sh, Sand water (0.1 μg L−1). Values of leached metabo(homogenous); Sp, Sand (planar x-beds); Sr, Sand (ripple x-laminated); St, Sand lites out of the root zone reached values above (trough x-beds); Sv, Sand (horizontal laminated).

Iversen et al.: Hydrogeological Relationships of Sandy Deposits

1915

advantage of using a distributed model for modeling leaching of pesticides on a larger scale, especially if the study is focusing on early warning of pesticide leaching. The need for developing indirect methods for up-scaling the lithofacies distribution from excavation scale to field scale or to an even larger scale is, therefore, present. A solution of this might be the use of ground-penetrating radars (GPR). For instance, Møller and Vosgerau (2006) performed a GPR survey on the same glaciofluvial deposits used in the current study. They concluded that radar stratigraphic analysis revealed realistic facies architecture models for the outwash plain and that GPR can be used for up-scaling lithofacies distributions in sedimentological settings.

Summary Eleven lithofacies from soils originating mainly from glaciofluvial deposits were described in terms of soil texture distribution, sorting, bedding style, and external Fig. 6. (a) and (c) Fitted van Genuchten/Mualem soil water characteristics and hydraulic conductivity boundaries. Facies ranging from gravelly functions for the investigated lithofacies and (b) and (d) the resulting hydrofacies. Values to sandy and diamicts were found. With for Hydrofacies 4 were found by including the four facies (Sh, Sp, St, and Sv) in the fitting respect to their soil hydraulic properties, procedure of the RETC code. Gd, Gravel (homogeneous or poorly layered); Sr, Sand (ripple the 11 different lithofacies formed four x-laminated); Dmh, Diamicts (matrix supported, homogeneous); Sh, Sand (homogenous); Sp, Sand (planar x-beds); St, Sand (trough x-beds); Sv, Sand (horizontal laminated). different hydrofacies having relatively homogenous properties. Differences in probably also related to the uncertainties in the modeling paramthe hydraulic properties were mainly related to differences in the eters related to the properties of the pesticide. Uncertainties in soil texture distribution and to a lesser extent to differences within model parameters with respect to pesticide leaching are known the depositional structure of the lithofacies. Two large outcrop to be high. For instance, Dubus and Brown (2002) simulated sections from one of the locations (a gravel pit) located near the leaching of two hypothetical pesticides and found that the coefterminal moraine of the former Weichsel glacier were used for the ficient of variation (CV) for pesticide losses ranged from 60 to two-dimensional modeling. Modeling results revealed that the 150%, which were in sharp contrast with percolation volumes spatial distribution of sedimentary bodies aff ected water flow and (CV around 6–7%). Finally, the process description in the apthe leaching of pesticides especially in the fi rst part of the simulaplied model may not be totally representative. This is related to tion period. In the last part of the simulation period a diminishing the model structure error (Refsgaard et al., 2006) or related to the eff ect was observed along the profi le. Th is was probably related representation of the modeled domain, which in the present artito the eff ect that an increased accumulated amount of leachate at cle is supported by observed data. If in our work a direct compar3-m depth applied over time resulted in a spatial averaging out of ison against measured field data had been possible, a calibration the diff erences in hydraulic properties of the hydrofacies as well as of the model would most likely have been necessary. Therefore, of small scale heterogeneities at that depth leading to a higher and it is not possible to compare the modeled values to the EU limit more evenly distributed concentration profile. This indicated that value for drinking water. However, recent work in Denmark on for the studied outcrop sections, differences in the hydraulic propidentifying areas vulnerable to pesticide leaching (Nygaard et al., erties were too small to create any major preferential flow patterns. 2004) points out that sediment types such as glaciofluvial deposits belong to the most vulnerable sandy areas in Denmark. Despite the similarity between the two hydrofacies (1 and 4), a response characteristic on the flow of water and leaching of metribuzin in the two outcrop sections could be observed from the model output after a simulation period of 3 yr. The response characteristic was reflected in the variation of the pesticide concentration differences along the 3-m line. This emphasizes the

Acknowledgments

1916

Journal of Environmental Quality • Volume 37 • September–October 2008

This work was supported by the KUPA project (concept for identifying areas where shallow aquifers are vulnerable to pesticide contamination), financed by the Danish Ministry of Environment and Energy. The authors wish to thank the technical staff at the Danish Inst. of Agricultural Sciences

(S.T. Rasmussen, M. Koppelgaard, and B.B. Christensen) for their valuable work in the field and in the laboratory.

References Alfnes, E., W. Kinzelbach, and P. Aagaard. 2004. Investigation of hydrogeologic processes in a dipping layer structure: I. The flow barrier effect. J. Contam. Hydrol. 69(3–4):157–172. Anderson, M.P. 1989. Hydrogeologic facies models to delineate large-scale spatial trends in glacial and glaciofluvial sediments. Geol. Soc. Am. Bull. 101(4):501–511. Anderson, M.P., J.S. Aiken, E.K. Webb, and D.M. Mickelson. 1999. Sedimentology and hydrogeology of two braided stream deposits. Sediment. Geol. 129(3–4):187–199. Bersezio, R., A. Bini, and M. Giudici. 1999. Effects of sedimentary heterogeneity on groundwater flow in a Quaternary pro-glacial delta environment: Joining facies analysis and numerical modelling. Sediment. Geol. 129(3–4):327–344. Christiansen, L. 2003. Eksperimentelle og numeriske undersøgelser af strømning og transport ved VAP lokaliteter (in Danish). M.S. thesis. Univ. of Copenhagen, Geological Inst., Copenhagen, Denmark. Dubus, I.G., and C.D. Brown. 2002. Sensitivity and first-step uncertainty analyses for the preferential flow model MACRO. J. Environ. Qual. 31:227–240. Heinz, J., S. Kleineidam, G. Teutsch, and T. Aigner. 2003. Heterogeneity patterns of Quaternary glaciofluvial gravel bodies (SW-Germany): Application to hydrogeology. Sediment. Geol. 158(1–2):1–23. Iversen, B.V., M. Koppelgaard, and O.H. Jacobsen. 2004. An automated system for measuring air permeability and hydraulic conductivity in the laboratory on large soil cores. DIAS report- Plant Production 111. Danish Inst. of Agricultural Sciences, Tjele, Denmark. Jensen, K.H., K. Bitsch, and P.L. Bjerg. 1993. Large-scale dispersion experiments in a sandy aquifer in Denmark–Observed tracer movements and numerical-analyses. Water Resour. Res. 29(3):673–696. Ju, S.H., and K.J.S. Kung. 1993. Finite element simulation of funnel flow and overall flow property induced by multiple soil layers. J. Environ. Qual. 22:432–442. Ju, S.H., and K.J.S. Kung. 1997. Impact of funnel flow on contaminant transport in sandy soils: Numerical simulation. Soil Sci. Soc. Am. J. 61:409–415. Kjær, J., P. Olsen, H.C. Barlebo, T. Henriksen, R.K. Juhler, F. Plauborg, R. Grant, P. Nygaard, and L. Gudmundsson. 2005a. The Danish Pesticide Leaching Assessment Programme, Monitoring results May 1999–June 2004. Geological Survey of Denmark and Greenland, Danish Inst. of Agricultural Sciences, National Environmental Research Inst., Copenhagen, Denmark. Kjær, J., P. Olsen, T. Henriksen, and M. Ullum. 2005b. Leaching of metribuzin metabolites and the associated contamination of a sandy Danish aquifer. Environ. Sci. Technol. 39(21):8374–8381. Klingbeil, R., S. Kleineidam, U. Asprion, T. Aigner, and G. Teutsch. 1999. Relating lithofacies to hydrofacies: Outcrop-based hydrogeological characterisation of Quaternary gravel deposits. Sediment. Geol. 129(3– 4):299–310. Klute, A. 1986. Water retention: Laboratory methods. p. 635–662. In A. Klute (ed.) Methods of soil analysis, Part 1. Physical and mineralogical methods. ASA and SSSA, Madison, WI. Klute, A., and C. Dirksen. 1986. Hydraulic conductivity and diffusivity: Laboratory methods. p. 687–734. In A. Klute (ed.) Methods of soil analysis, Part 1. Physical and mineralogical methods. ASA and SSSA, Madison, WI. Kung, K.J.S. 1990. Preferential flow in a sandy vadose zone: II. Mechanism and implications. Geoderma 46(1–3):59–71. Madsen, H.B., A.H. Nørr, and K.Aa. Holst. 1992. The Danish Soil

Iversen et al.: Hydrogeological Relationships of Sandy Deposits

Classification. Atlas of Denmark. I,3. The Royal Danish Geographical Society, Copenhagen, Denmark. Miall, A.D. 1977. A review of the braided-river depositional environment. Earth-Sci. Rev. 13(1):1–62. Miall, A.D. 1983. Glaciofluvial transport and deposition. p. 168–183. In N. Eyles (ed.) Glacial geology. Pergamon Press Ltd., Oxford, UK. Møller, I., and H. Vosgerau. 2006. Testing ground-penetrating radar for resolving facies architecture changes–A radar stratigraphic and sedimentological analysis along a 30 km profile on the Karup Outwash Plain, Denmark. Near Surf. Geophys. 4(1):57–68. Mualem, Y. 1976. A new model for predicting the hydralic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522. Nygaard, E., V. Ernstsen, C.S. Jacobsen, O.H. Jacobsen, R. Juhler, P. van der Keur, S.E. Olesen, J. Rasmussen, P. Rosenberg, and H. Vosgerau. 2004. Pesticide leaching in Danish groundwater: Identification of vulnerable areas. Geological Survey of Denmark and Greenland Bull. 4:25–28. Olesen, J.E., and T. Heidmann. 1990. EVACROP. Et program til beregning af aktuel fordampning og afstrømning fra rodzonen, version 1.00 (In Danish). Arbejdsnotat nr. 9. Dep. of Agrometeorology, Danish Inst. of Plant and Soil Science, Tjele, Denmark. Refsgaard, J.C., J.P. van der Sluijs, J. Brown, and P. van der Keur. 2006. A framework for dealing with uncertainty due to model structure error. Adv. Water Resour. 29(11):1586–1597. Simunek, J., M. Sejna, and M.Th. van Genuchten. 1999. The HYDRUS2D software package for simulating two-dimensional movement of water, heat, and multiple solutes in variably saturated media. Version 2.0, IGWMC-TPS-53. International Ground Water Modeling Center, Colorado School of Mines, Golden, CO. Torp, S. (ed.). 2005. Ekstramarginale smeltevandsaflejringer af Sen Weichsel alder inden for alluvialkeglen: Basisdata fra Karup og Tinglev hedesletter (In Danish). Rep. No. 4. Geological Survey of Denmark and Greenland, Danish Inst. of Agricultural Sciences, Copenhagen, Denmark. van den Elsen, E., J. Stolte, and G. Veerman. 1999. Three automated laboratory systems for determining the hydraulic properties of soils. p. 329–340. In M.T. van Genuchten, F.J. Leij, and L. Wu (ed.) Characterization and measurement of the hydraulic properties of unsaturated porous media, Part 1. U.S. Salinity Lab., USDA-ARS. Dep. of Environmental Sciences, Univ. of California, Riverside, CA. van der Keur, P., M.P. Hag, P. Rosenberg, J.J. Rasmussen, B.V. Iversen, and C.D. Børgesen. 2004. Identification of vulnerable sandy areas with respect to pesticide leaching using a pedotransfer function-process based model approach. p. 349–355. In P. Agaard et al. (ed.) Proceedings of the International Workshop, Rome, 5–7 May 2004. COST Action 629: Fate, impact, and indicators of water pollution in natural porous media. La Goliardica Pavese, Pavia, Italy. van Genuchten, M.T. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898. van Genuchten, M.T., F.J. Leij, and S.R. Yates. 1991. The RETC code for quantifying the hydraulic functions of unsaturated soils Version 1.0. EPA Rep. 600/2-91/065, U.S. Salinity Lab., USDA-ARS, Riverside, CA. Vereecken, H., U. Döring, D.J. Kim, J. Feyen, C. Mouvet, C. Moreau, P. Burauel, and K. Hucker. 1995. The fate and mobility of pesticides at laboratory, lysimeter, and field scale: A European project. p. 147–154. In A. Walker et al. (ed.) Proceeding of a symposium held at the Univ. of Warwick, Coventry, UK, 3–5 Apr. 1995. Pesticide movement to water. Monogr. No. 62. British Crop Protection Council, Farnham, UK. Walter, M.T., J.-S. Kim, T.S. Steenhuis, J.-Y. Parlange, A. Heilig, R.D. Braddock, J.S. Selker, and J. Boll. 2000. Funneled flow mechanisms in a sloping layered soil: Laboratory investigation. Water Resour. Res. 36(4):841–849.

1917