Modeling subsurface transport in extensive glaciofluvial and littoral sediments to remediate a municipal drinking water aquifer

Few studies have been carried out that cover the entire transport process of pesticides, from application at the soil surface, through subsurface transport, to contamination of drinking water in esker aquifers. In formerly glaciated regions, such as Scandinavia, many of the most important groundwater resources are situated in glaciofluvial eskers. The purpose of the present study was to model and identify significant processes that govern subsurface transport of pesticides in extensive glaciofluvial and littoral sediments. To simulate the transport processes, we coupled a vadose zone model at soil profile scale to a regional groundwater flow model. The model was applied to a municipal drinking-water aquifer, contaminated with the pesticide-metabolite BAM (2,6-dichlorobenzoamide). At regional scale, with the combination of a ten-meter-deep vadose zone and coarse texture, the observed concentrations could be described by the model without assuming preferential flow. A sensitivity analysis revealed that hydraulic conductivity in the aquifer and infiltration rate accounted for almost half of the model uncertainty. The calibrated model was applied to optimize the location of extraction wells for remediation, which were used to validate the predictive modeling. Running a worst-case scenario, the model showed that the establishment of two remediation wells would clean the aquifer in four years, compared to nine years without them. Further development of the model would require additional field measurements in order to improve the description of macrodispersion in deep, sandy vadose zones. We also suggest that future research should focus on characterization of the variability of hydraulic conductivity and its effect on contaminant transport in eskers. Correspondence to: M. Bergvall (martin.bergvall@tyrens.se)


Introduction
An increasing number of groundwater aquifers and drinkingwater supplies across the world are being contaminated by chemicals from anthropogenic sources (Erickson, 2002;EEA, 1999;USEPA, 1999).Pesticides and their degradation products are some of the most common contaminants reported in European countries (Lapworth et al., 2006;Cerejeira et al., 2003;Papastergiou and Papadopoulou-Mourkidou, 2001;Schiavon et al., 1995) and in the US (Barbash et al., 2001;Kolpin et al., 1998).In Sweden, there are about 1900 municipal drinking-water plants that draw their supplies from natural or artificial groundwater.As part of a national monitoring program, 28 % of the groundwater supplies investigated between 2005 and 2008 were identified as being contaminated with pesticides (CKB, 2009;Törnquist et al., 2007).BAM (2,6-dichlorobenzoamide) was found to be the most common contaminant, reported in 60 % of the contaminated groundwater supplies.Similar levels of contamination have been found in the Netherlands, Germany, and France (Morvan et al., 2006;Versteegh and te Biesebeek, 2003;Wolter et al., 2001).In Denmark, for example, BAM has been detected in 21 % of all groundwater wells investigated (Clausen et al., 2004) and approximately one hundred groundwater extraction wells were recently closed due to BAM contamination (Sørensen et al., 2009).
BAM is a hydrophilic metabolite of dichlobenil (2,6dichlorobenzonitrile) and is highly persistent in soil (Cox, 1997).Dichlobenil has largely been used as herbicide to control various types of weeds.The use of the product Totex, which contains the active substances dichlobenil and atrazine, has been banned in Sweden since 1991, and in many other countries (KEMI, 2008).The toxicity of dichlobenil, which is classified as slightly to moderately toxic, has been investigated by several authors (Cox, 1997;Linders et al., 1994).However, relatively little is known about the toxicity and health effects of its degradation product BAM, although the available data indicate that its potential carcinogenicity does not exceed that of dichlobenil (USEPA, 1998).The potential effects of long-term exposure to low concentrations of either compound have not been thoroughly investigated.
The large number of contaminated groundwater-supplies highlights the need for tools to evaluate the subsurface processes controlling contaminant transport.The hydrogeology of groundwater aquifers is in general complex, and distributed numerical modeling is often required.Depending on the extent to which a model can be validated using field data, it can be used to predict the behavior of contaminants over reasonably long time-scales and under various conditions (Hassan, 2004).Numerical models can also be used to predict the potential of remediation strategies at a specific site and aid decision making (Zhou et al., 2004;Prommer et al., 2000).
Most numerical models are optimized to solve contaminant transport problems in either the vadose or groundwater zone, however, models have been proposed which couple transport in both zones (Weill et al., 2009;Panday and Huyakorn, 2008;Twarakavi et al., 2008;Trefry and Muffels, 2007).For pesticide transport, one-dimensional vadose zone models and three-dimensional groundwater models have been coupled and applied at local (Stenemo et al., 2005;Jorgensen et al., 2004) and regional scales (Herbst et al., 2005;Christiansen et al., 2004).These models were developed for sites with fine-textured soils and shallow groundwater, where macropore flow is a significant transport process.
An esker is a narrow ridge of gravelly and sandy glacial material originally deposited by a stream in an ice-walled valley.In esker aquifers, which are generally extensive in depth and range few regional-scale studies have considered the whole transport process, from application of pesticide at the soil surface, through transport in the vadose zone, to contamination of drinking water (Köhne et al., 2009).The lack of models for glaciofluvial sediments at this scale is problematic since many of the most important groundwater resources in recently glaciated regions of the world are situated in such aquifers.
The objectives of this study were: (1) to model and identify significant processes that govern subsurface transport in extensive glaciofluvial and littoral sediments, focusing on a BAM-contaminated esker aquifer, which supplies a municipal drinking-water utility; and (2) to use knowledge of these processes to design and test a remediation strategy and validate the results from the model.

Modeling approach
Simulation of pesticide movement was performed in two steps.First, transport in the vadose zone was modeled at soil profile scale with HYDRUS-1D ( Šimunek et al., 1998).HYDRUS-1D is a 1-D finite-element model which accounts for a variety of processes, including soil-water flow, solute transport, heat transport, sorption, degradation, volatilization, crop uptake and surface runoff.To simulate transport in groundwater at the regional scale, the calculated pesticide concentration from the bottom of the vadose zone model was used as the input, through the recharge parameter, to the 3-D finite-difference model MODFLOW (Harbaugh et al., 2000) coupled with MT3DMS (Zheng and Wang, 1999).MOD-FLOW calculates flow and MT3DMS, which is designed to interface with the MODFLOW code, performs transport calculations (Fig. 1).
HYDRUS-1D numerically solves the Richards equations for variably saturated water-flow and advection-dispersion type equations for solute transport.The water-flow equation solved by HYDRUS-1D (for flow in the vertical direction) is: where ψ is the pressure head; θ is the volumetric water content; t is time; z is the spatial coordinate (positive upward); q s is the sink term; and K is the unsaturated hydraulic conductivity.The van Genuchten-Mualem model was used to calculate soil hydraulic properties (van Genuchten, 1980).The solute transport equation with first-order reaction rate solved by HYDRUS-1D is: where s and c are the solute concentrations in the solid and liquid phases, respectively; q is the volumetric flux density; ρ is the soil bulk density; θ is the volumetric water content; t is time; z is the spatial coordinate (positive upward); D is the dispersion coefficient for the liquid phase; and µ is a first-order rate constant.
The transient flow equation solved by MODFLOW is: where h is the hydraulic head; K ii is the tensor of hydraulic conductivity; x i is the distance along the respective Cartesian coordinate axis; q s is the sink/source term; S S is the specific storage; and t is time.The advection-dispersion equation solved by MT3DMS for one species is:

∂ (p
where C is the solute concentration; p k is the kinematic porosity; t is time; x i is the distance along the respective Cartesian coordinate axis; D ij is the tensor of hydrodynamic dispersion; v i is the seepage velocity; q s is the sink/source term related to the aquifer; C s is the concentration of the sink/source flux; ρ b is the bulk density; S is the sorbed concentration; and R n is the chemical reaction term. To evaluate the robustness of the model's predictions, a sensitivity analysis was undertaken with respect to the physical and chemical parameters included in the model (Eqs. 1 up to 4).Considered parameters were the van Genuchten-Mualem parameters, infiltration rate, thickness of the vadose zone, dispersivity, sorption and degradation of the contaminants, and hydraulic conductivity in the aquifer.To calculate the effect that a parameter has on the model outcome, the parameters in the calibrated models were perturbed with changes of ±2 to 10 %.The equation used to compare the perturbed models was: where RMS is the root-mean-square deviation in the prediction; c cali is the solute concentration at the target in the calibrated model at time step i; c peri is the solute concentration at the target in the perturbed model at time step i; and N is the number of time steps.

Field site description
The city of Umeå is located in northern Sweden (lat.63 • 83 N, long.20 • 26 E) and has a population of 115 000.Groundwater is an important source of municipal drinking-water for this city.Unfortunately, the concentration of BAM in the groundwater exceeds the European drinkingwater standard for pesticides, 0.1 µg l −1 (EU, 2006).Groundwater measurements began in 1998 and, since then, BAM has been at concentrations exceeding the standard, by up to 22 times, in two extraction wells.A former forest-plant nursery situated in Piparböle, about 1.5 km up-gradient from the water-extraction wells (Fig. 2), is the source of the contaminant (Bergvall et al., 2007).The nursery, which was active between 1956 and 2002, cultivated mainly spruce and pine seedlings.Several kinds of pesticides were used at the nursery, but only dichlobenil and its degradation product BAM were considered in this study since there is great concern over BAM contamination owing to its high mobility.Records from the nursery suggest that Totex, with the active substances dichlobenil and atrazine, was initially applied during the transition from growing bareroot seedlings on open land to their production in small containers on open land.The transition between production methods began during the 1970s and records suggest that the first application of Totex was in 1976.Totex was applied in granules, mainly along the plant beds and at the short sides where the irrigation system was located.According to historical aerial photographs, this area was estimated to cover 4 % of the total 120 000 m 2 nursery area.The recommended application rate of the active substance was 0.16 to 1.7 g m −2 .
The former nursery is situated on well-sorted, unstructured, littoral sand and lenses of silt.The soil has not been developed but mechanically cultivated, and the organic carbon content is low.The littoral sand is in contact with a north-south-oriented esker, which consists of glaciofluvial sand and gravel containing large quantities of groundwater.The general direction of groundwater movement is southwards, although local-scale variation in this movement does occur (Fig. 2).The groundwater table at the nursery is at a depth of approximately 10 m, but at the extraction wells its depth is 23 to 27 m, depending on the pumping rate.
The esker is partly covered by thin layers of silt and clay, especially along the sides.Layers of silt and clay also underlie the northern and southern lakes (Fig. 2), which explains why the groundwater head at the southern end of the southern lake is found 27 to 28 m below the surface water head.However, the esker is also partly located beneath the silty and clayey beds of the lakes.

Data sampling -unsaturated zone
In 2003, the soil surface at the former nursery was screened for dichlobenil and BAM.At one identified hotspot (Fig. 2), soil samples were collected at 40 cm intervals, down to the groundwater table, in November 2003.Details of the sampling have been reported previously (Bergvall et al., 2007) and the results are summarized in  method (105 • C for 12 h), and then measured on thirteen occasions with a depth moisture probe (Nucletronics ApS, Borre, Denmark) (Fig. 3).
In February 2009, a total of thirty soil samples were collected.They originated from five locations and were taken every 50 cm, down to a depth of 5 to 7 m (Fig. 2).The samples were frozen (at −18 • C) until analysis.Soil samples were analyzed with liquid chromatography-tandem mass spectrometry (LC-MS/MS, Micromass Quattro Ultima coupled to an Agilent 1100 HPLC).The concentration of total organic carbon (TOC) was calculated as the percentage weight lost on ignition, in accordance with the SS-EN 12879 method (SIS, 2000).The soil profile that was sampled in 2003 and the five soil profiles that were sampled in 2009 showed similar layering of sandy sediments.

Data sampling -saturated zone
The hydraulic groundwater head was monitored in 82 observation wells (Fig. 2) in October 2005 and August 2007.The monitoring took one week in each case.Fifty five of the observation wells were also used for sampling (Fig. 2).The groundwater collected was subsequently analyzed for BAM concentration using LC-MS/MS.The maximum number of samples analyzed from a single observation well at the nursery was 46 (Fig. 4), whereas in peripheral wells only one sample was analyzed (Table 2, Fig. 5).In total, 101 samples have been analyzed from the two extraction wells since 1998 (Fig. 6).One problem with the observations shown in Figs. 4 and 6 is that many different people from the local authority collected these samples.This means there is uncertainty as to whether the wells were purged adequately prior to sampling.Therefore, Fig. 4 also shows a moving average, which is the mean value of the closest four observations, before and after each observation.This means that the integral of the observed concentrations should over time equal the integral of the moving average concentrations.In addition, the extraction wells have not operated on a regular basis, which may have affected the sampling mixture of the collected water (Fig. 6).Since 2007, the wells have been running continuously for remediation purposes and the sampling frequency has been intensified.
The stream that drains from the southern lake crosses the aquifer above the groundwater surface.We estimated the drainage from the stream to the aquifer by measuring the stream discharge above and below the aquifer-crossing, Table 1.Analysis of soil samples taken in six, 5 to 10 m deep, soil profiles from the nursery site, including data from Bergvall et al. (2007), and slug tests carried out in the aquifer (Fig. 2).Except for the hydraulic conductivity (K gw ) measurements from slug tests, all data refer to the vadose zone.For each parameter the mean, standard deviation (SD), observed range and number of observations are reported.When the concentration was less than the detection limit (0.1 µg kg −1 ) 0.05 µg kg −1 was used in the calculation of the mean and standard deviation.using two water capacitance probes (TruTrack loggers WT-HR 64 K, New Zealand).The water-level probes were calibrated on four occasions in 2005, using an automated saltdilution method (Q-trace, Logotronics, Austria).On the first occasion, salt injections were also used at four other downstream locations (Fig. 2).The mean leakage was estimated to be 11 l s −1 , which is 7 % of the mean discharge rate (SMHI, 2009).Two pumping tests were conducted in the study area.One extraction well, situated 1100 m north of the model boundary (Fig. 2), was pumped for 127 days at 112 l s −1 .The two BAM-contaminated extraction wells, located 60 m apart, were pumped for 10 days at 83 l s −1 .Both tests were performed in unconfined parts of the esker aquifer.Results from the 127-day test indicated the hydraulic conductivity was between 90 and 290 m d −1 .The specific yield was between 7.4 and 9.1 %.The 10-day test indicated a hydraulic conductivity of 90 m d −1 and a specific yield between 7.0 and 20.1 %.
Hydraulic conductivity was also determined using slug tests in the observation wells (Fig. 2) (Hvorslev, 1951).The tests were carried out in a wide range of sediments including esker material, sand, silt and till (Table 1).
At the nursery, a local-scale, salt-injection study was conducted between August and October 2003.Salt water was injected into the groundwater in one well and transported by natural gradient flow to another well 2 m away.At a rate of 0.38 to 0.98 l min −1 1200 l of salt water, with an electrical conductivity of 147 m S m −1 , was injected.The electrical conductivity was logged every 30 min at three levels in both wells (at depths of 0.1 m, 1.4 m and 1.9 m below the groundwater table ).Assuming Gaussian distribution in the flow direction, the longitudinal dispersion was calculated by an analytical solution of the Fickian advection-dispersion model of solute transport, and the resulting values ranged between 0.2 and 1.5 m, at the 2 m scale.

Vadose zone model
The vadose zone model was only applied to the nursery area (Fig. 2).Since there was little variation in groundwater level within the nursery area (standard deviation = 0.06 m, Table 1), we assumed no feedback from the groundwater model to the vadose zone model.A 1-D model was motivated since there was no evidence of tilting layers under the soil surface, which can cause lateral flow.The slope of the surface was small (0.8 %) and the soil-profile model (Fig. 1), which was assumed to be representative of the total nursery Hydrol.Earth Syst.Sci., 15,[2229][2230][2231][2232][2233][2234][2235][2236][2237][2238][2239][2240][2241][2242][2243][2244]2011 www.hydrol-earth-syst-sci.net/15/2229/2011/ Table 2. Concentrations of BAM for the four zones (Fig. 5) and three time intervals.For each zone and time interval, the mean, standard deviation (SD), range and number of samples are reported.When the concentration was less than the detection limit (0.01 µg l −1 ) 0.005 µg l −1 was used in the calculation of the mean and standard deviation.area, covered a depth of 10.2 m, from the soil surface to the groundwater table.It was assumed that there was no surface runoff because of the flatness and texture of the soil surface, which was classified as medium to coarse sand.
The details of the vadose zone model (Larsson and Grip, 2004) were as follows.The conceptual model included variable texture and water content, at different depths relative to the groundwater table.A finite-element grid was designed for modeling flow and transport.The grid was narrowly spaced at 2 cm in the vertical direction, which is sufficient resolution for valid calculations of sandy soils (Vogel and Ippisch, 2008).Ten soils with various properties were defined to describe the vadose zone, ranging from fine sandy silt to coarse sandy medium sand (Table 3).
Degradation of dichlobenil and BAM was described by first order kinetics.Conditioned on the total mass balance, the half-lives were calibrated to the values given in Table 3.The values agree with those of a Danish study (Clausen et al., 2007), which also demonstrated that various redox conditions do not affect degradation of dichlobenil.As in other studies, they found strong evidence that the degradation to BAM is a microbially catalyzed process.Microorganisms are most likely to be found in the upper layers of the soil and our shortest half-lives were set close to the soil surface.Clausen et al. (2007) could not correlate degradation of dichlobenil to BAM with sorption, water chemistry or composition of soils or sediments.
It was assumed that the Langmuir sorption isotherm was applicable.However, as the linear Freundlich isotherm produced similar results at depths greater than 0.9 m, the Langmuir isotherm was not used at these depths in order to save computing time.Clausen et al. (2004) found a linear correlation between the total organic carbon content (TOC) and sorption of dichlobenil and BAM as long as the TOC > 1.0 g kg −1 : where K d dich and K d BAM are the sorption coefficients of dichlobenil and BAM, respectively; and f oc is the weight fraction of TOC.The sorption coefficients for dichlobenil and BAM were initially calculated in accordance with Eqs. ( 6) and ( 7), but were later calibrated to improve the correlation with measured concentrations of dichlobenil in the solid and dissolved phases (Table 3).About 20 cm of the topsoil had been excavated before we sampled the soil, and since there had been vegetation we assumed about 3 to 5 times higher TOC-content in the topsoil layer (Table 3) compared to our analyzed maximum TOC (Table 1).
We assumed high macrodispersivity values to account for possible forms of preferential flow.Dispersivity, α vz , was initially assumed to be 10 % of the soil-profile depth, but was calibrated to 5 m except for the top meter.The dimensionless Henry coefficient was set to 3.5 × 10 −4 (Linders et al., 1994).
Since the pesticide was very expensive, we assumed that only 0.17 g m −2 was applied, twice a year between 1976 and 1983, and once a year between 1984 and 1991.This means that 19 kg of the active substance, dichlobenil, was applied to 4 % of the total nursery area between 1976 and 1991; this is considered to be a conservative estimate (Bergvall et al., 2007).However, since the soil has been repeatedly mixed over the years, we assumed the same load of pesticide over the entire nursery area (i.e. 4 % of 0.17 g m −2 ).The vadose zone model was calibrated with data sampled at the identified hotspot (Fig. 2).Considering the higher dichlobenil concentrations found at the hotspot, the application rate was assumed to be 1.7 g m −2 in the calibration model.
The simulation started on 1 September 1975 and the calibration model was run until 1 November 2003.The prediction model was run for 45 years until 31 December 2020.Daily climate data, including precipitation, temperature, wind velocity, relative humidity and global radiation, were collected from a weather station located 10 km southeast of the nursery (SMHI, 2009).A temperature-index approach was used to model snow cover and snow melt (Anderson, 1976), and potential evaporation was calculated according to the Penman equation (Penman, 1948).To generate input data for the predictive modeling, the climate data from the preceding eight years were used; for example, climate data from 2003 were used for 2011.During the simulation time, the maximum calculated infiltration rate was 40 mm da −1 , and occurred a day when the precipitation was 42 mm d −1 .Since this is 10 3 times smaller compared to the saturated hydraulic conductivity at the soil surface (Table 3), and rain with an intensity of 20 mm per half an hour occurs only once in ten years at the field site (SMHI, 2009), daily time resolution of input data was assumed to be satisfactory.Consequently, we considered matrix flow to be the dominant process.
The water flow and the solute transport models were calibrated sequentially.The van Genuchten-Mualem parameters including saturated hydraulic conductivity were calibrated against water contents analyzed by the oven-drying method, and then validated with the moisture probe data obtained at the hotspot (Fig. 3).The mean error of the total water balance was 1.3 %.The root-mean of the squared differences between the calculated and observed water contents obtained by the oven-drying method was normalized by the observed range, giving a calibration error of 10 %.The root-mean of squared differences for the validation data gave a corresponding error of 24 %.The infiltration rate and the generated bottom water flux at 10.2 m depth are shown in Fig. 7.The figure shows a seasonal variation of the bottom flux.
The solute component of the model was calibrated to achieve the best match between the calculated and observed hotspot-concentrations of dichlobenil in the solid and solute phases (Table 1, Fig. 8).The sensitivity analysis that was undertaken, Eq. ( 5), helped to identify the important parameters to be fitted.Since the observed concentrations had high variability the obtained fit to simulated concentrations was only in the same order of magnitude.The root-meansquares of the residuals of the adsorbed and liquid phases were normalized by the observed concentration range giving calibration errors of 38 % and 39 %, respectively.Concentrations and mass fluxes of BAM that were generated by the prediction model at the bottom of the vadose zone are shown in Fig. 9.
At the bottom of the vadose zone the mass fluxes of BAM corresponded with the variation of water fluxes (Fig. 7), and resulted in stable BAM concentrations (Fig. 9).In comparison with the mass fluxes of BAM from the calibration model, the prediction model generated correspondingly lower mass fluxes equivalent to the lower application rate, i.e. in the predictions it was valid to use the calibration model with a lower application rate.

Groundwater model
A 3-D, block-centered, finite-difference grid was designed to model flow and transport in the saturated zone.In the horizontal direction, a total of 73 columns and 153 rows (uniformly spaced at 20 m) were used for flow simulation, covering an area of 1.45 km by 3.05 km.The aquifer was divided into two layers, with the upper layer approximately 2 m below the groundwater table.
Boundaries were identified as general heads, specified heads, no-flows and rivers (Fig. 10).The conductances of the general head boundaries were calibrated to give inflows corresponding to the net infiltration rate in the northeast and to flow estimates from the two pumping tests (conducted north of the northwestern boundary and at the extraction wells); the inflow rate was calculated as 83 l s −1 in the northwest and 0.4 l s −1 in the northeast.The conductance of the southern The assumed application rate of dichlobenil was 0.17 g m −2 , which was ten times lower than the rate in the calibrated model.The BAM concentrations from the vadose zone model were used as input to the groundwater model.stream, where the stream flows above the esker, was calibrated to match stream leakage.The conductance of the eastern river boundary was set to permit small inflows and outflows to the aquifer.
In accordance with records from the municipality, the extraction rate of the two contaminated wells varied between 0 and 61 l s −1 .Four different recharge rates were defined for the model: most areas had a constant net infiltration rate of 300 mm d −1 .The constant rate at the nursery field was justified by the 10 m thick vadose zone, which smoothes out variations in the net infiltration rate and changes in the groundwater level (Fig. 7, Table 1).The BAM concentrations generated by the vadose zone model had also small variations from day to day (Fig. 9).The lakes were assumed to give no recharge to groundwater (Fig. 2).Areas in the vicinity of the lakes, mainly west of the southern lake, consist of two aquifers, where the upper one (Fig. 2, consisting of silt and clay) has a leakage to the lower one (consisting mainly of sand).The upper aquifer has primarily an outflow to the lakes and is not connected to the contaminated area.A leakage from the upper aquifer was, therefore, modeled as recharge to the lower aquifer, assuming a net infiltration rate of 30 mm d −1 .The erosion-cuts along the streams (Fig. 2) consist mainly of sand and the recharge rate was assumed to be 15 mm d −1 .
Initially, hydraulic conductivity was quantified in accordance with the results from the pumping and slug tests.The extent of the zones of hydraulic conductivity were derived from borehole data, based on 91 borings (Wikner and Öhlén, 1994), and ground-penetrating radar measurements (RA-MAC/GPR, 100 MHz antenna, Malå Geoscience, Sweden), using prior knowledge of how glaciofluvial sediments have been geologically deposited (Fig. 10).
The bedrock surface defines the bottom of the model (Fig. 10), except in small areas where silt or clay is deposited beneath more coarse-grained textures.The calibration of hydraulic conductivity zones in the two layers was conducted in two steps to match the calculated and observed heads at the observation wells.First, assuming no pumping at the extraction wells, the model was calibrated to the observations from 2005.Secondly, with the total pumping rate set to 44 l s −1 in the extraction wells, the model was calibrated to the heads observed in 2007.The calibrated flow model had a mean head error of 0.11 m.The calibration error, defined as the root-mean of the squared head differences normalized by the observed head range, was 1.6 %, which is less than the calibration target of 5 %.
A transient validation of the model was conducted using data from the 10-day pumping test performed in the extraction wells.The validation was also used to confirm that the mean of the specific yield evaluated from the pumping test was accurate.Since the aquifer has coarse texture and is unconfined at the location where the pumping test was conducted, the kinematic porosity was estimated from the mean specific yield (0.13).
According to Clausen et al. (2007), only insignificant degradation of dichlobenil and BAM has been observed in groundwater.Since dichlobenil has not been detected in our aquifer, only transport of BAM was modeled.The solute concentration data from the HYDRUS-1D model was evenly distributed in the grid cells corresponding to the open land area of the former nursery (Fig. 2).Time was divided into half-year intervals, which provided sufficient resolution to respond to changes in the input of the contaminant and mass conservation.The simulation started on 1 January 1976 and the calibration model was run until 1 July 2009.The prediction model was run until 31 December 2020.
Calibration of dispersivity in the three dimensions was conducted to correlate the calibrated with the observed concentrations of BAM in the observation and extraction wells (Table 2).By considering the integral of the observed concentrations, the aim was to fit the concentrations in observation well #7705 (Fig. 2) and the extraction wells with the least-squares method.The root-mean-squares of the residuals of the observed and modeled BAM concentrations at the observation well and at the extraction wells, were normalized by the observed range of concentrations giving calibration errors of 15 % and 14 %, respectively.
Calibration of the transport model gave a longitudinal dispersivity of 5 m, a transverse horizontal dispersivity of 0.1 m and a vertical dispersivity of 0.04 m.This is consistent with the results from the local-scale, salt-injection study (α L between 0.2 and 1.5 m), and can be put into perspective with the empirical relationship between longitudinal dispersivity (α L ) and field scale (L), suggested by Schulze-Makuch ( 2005): Equation ( 8) gives a longitudinal dispersivity of 32 m at the field scale (L = 1500 m), 0.96 m at the grid-cell scale (L = 20 m), and 0.15 m at the scale of the salt-injection study (L = 2 m).Both the calibrated and the measured longitudinal dispersivity values were comparatively higher than the values obtained from Eq. ( 8).The relatively low transverse horizontal dispersivity and vertical dispersivity values are realistic since advection dominates the mass transport.

Remediation strategy and model validation
Since May 2007, water from the two contaminated extraction wells has been pumped into a nearby river.Because these wells are situated 1.5 km from the contamination source, this has not been a cost-effective way to remediate the area.
The contaminated vadose zone amounts to approximately 1.2 million m 3 of sand.Therefore, application of a treatment system to the vadose zone is unfeasible.
The calibrated model was used to design a more costeffective remediation strategy.Locations, dimension and the number of new extraction wells for remediation purposes were optimized.Two new extraction wells were established (Fig. 2) and a 7-day pumping test was conducted in order to validate the calibrated hydraulic conductivity values and predicted BAM concentrations.

Results
Considering the model was run at daily time resolution, the saturated hydraulic conductivity of the aquifer was found to be the most sensitive parameter for determining the predicted concentration of BAM in the extraction wells (Table 4).
Together with the infiltration rate these parameters account for almost half the model uncertainty.The degradation rate of dichlobenil, which was expected to be one of the most uncertain parameters, was found to be the sixth most sensitive parameter.Since the vadose zone model was applied on a coarse-grained unstructured soil with high saturated hydraulic conductivities exceeding normally occurring infiltration rates by orders of magnitude, the validity of the results from the sensitivity analysis might be extended to time scales smaller than one day.
The simulated and observed concentrations, at one observation well at the nursery and at the extraction wells, are shown in Figs. 4 and 6.Calibration to the moving average gave a correlation coefficient of 0.99, compared to 0.63 for the observed values (Fig. 4).The correlation coefficient between the observed and the simulated concentrations are 0.40 for the western extraction well (Fig. 6).The simulation indicated that the guideline limit was exceeded in 1985 and that the concentration reached a maximum at the end of 1995.
The simulated BAM plumes in groundwater from 1985, 1995 and 2012 are shown in Fig. 11.This figure shows how the plume moves in the direction of the extraction wells and that the plume gets progressively smaller towards the front.The main leakage takes place at the southern end of the treenursery area, but there is also a small amount of leakage in the northwest.The two plumes enter the esker and become a single plume before reaching the extraction wells.This scenario is supported by the BAM samples in the area (Table 2, Fig. 5).A worst-case scenario was run, based on the uncertainty of the observed and simulated concentrations at the extraction wells.The standard deviations of the residuals of observed and simulated concentrations were used to calculate a 98 % confidence limit (±0.095 µg l −1 ).In the worst-case scenario, it was assumed that the model underestimates the true concentration, corresponding to the 98 % confidence limit.With this condition it was predicted that the BAM concentration will fall below the guideline limit in 2019.
In the modeled scenario, 19 kg of the active substance, dichlobenil, was applied between 1976 and 1991.At the end of the predictive simulation in 2020, 0.8 kg remained in the vadose zone as either dichlobenil or BAM, 2.6 kg had left the vadose zone as dissolved dichlobenil or through volatilization of dichlobenil or degradation of BAM, and 15.6 kg had degraded from dichlobenil to BAM and entered the groundwater.Of the 15.6 kg that had entered the groundwater, 8.2 kg was predicted to have been extracted in the wells and 2.5 kg was predicted to have passed the extraction wells because they have been periodically closed.
To calculate transit times in the vadose zone model, the medians of the water content at each calculation node was used.The medians were received from daily output during one simulation year (1 July 2003 to 30 June 2004).At each calculation node the retardation factor and transit times of water and BAM could then be calculated, and resulted in a median transit time to groundwater of 4.3 years for water and 7.2 years for BAM.This corresponds to a median retardation factor of 1.7, although this value varies depending on moisture content, depth, and sorption coefficients at different depths.In order to calculate transit times with the calibrated groundwater flow model, particle tracking was conducted from the groundwater table of the nursery field.It gave a median transit time between the nursery field and the extraction wells of 2.7 years for water and 4.7 years for BAM, which also corresponds to a retardation factor of 1.7.The transit time of groundwater from the northern to the southern end of the nursery was calculated to be more than one year, implying that local variability in mass fluxes will be mixed at regional scale.
The results of the 7-day pumping test in the two newly established remediation wells (Fig. 2) were compared to the model predictions in Table 5.During the pumping test, groundwater was analyzed twice and the mean concentrations were −1 % and +33 %, compared to the model results.According to the calibrated model, the recharge rate of the nursery's catchment area is 2.6 l s −1 .During the test, the pumping rates in the two wells were 0.82 l s −1 and 0.27 l s −1 .Since the evaluation of the pumping test indicated that more groundwater could be extracted, the simulation was run with pumping rates of 1.5 l s −1 and 0.75 l s −1 .With the remediation wells running from July 2010, and imposing the uncertainty of observed and simulated concentrations at the contaminated extraction wells (98 % confidence limit ±0.095 µg l −1 ), it was predicted that the concentration of BAM in the extraction wells would decrease to less than half the EU guideline limit at the latest of 2014.This corresponds to a reduction of five years compared to the situation without the remediation wells.

Discussion
The results were received from a glaciofluvial esker aquifer contaminated by the herbicide degradation product BAM.Despite the herbicide has been modestly used, BAM has contaminated an aquifer holding large groundwater quantities, and particularly two extraction wells located about 1.5 km down-gradient from the source have concentrations that exceed the European drinking water standard for pesticides (0.1 µg l −1 ).At regional scale, the model was able to explain observed BAM concentrations in the groundwater, and to evaluate the effect of remediation strategies.
On the groundwater table underlying the source the breakthrough curve of BAM became wide (Fig. 4).Similar results have been obtained in previous studies on pesticide leaching in unstructured soils (Kjaer et al., 2005;Brown et al., 2000).Since sandy soils naturally contain a high number of large pores, they permit infiltration when exposed to heavy rainfall.A lack of macropores implies that matrix flow is the dominant process, which causes a wide breakthrough curve.This helps to explain why the BAM contaminant was predicted at levels exceeding the drinking water standard at the two extraction wells for three to four decades.
Hydrol.Earth Syst.Sci., 15,[2229][2230][2231][2232][2233][2234][2235][2236][2237][2238][2239][2240][2241][2242][2243][2244]2011 www.hydrol-earth-syst-sci.net/15/2229/2011/In groundwater the results showed a limited lateral dispersal of the contaminant BAM, as it moves from the source in littoral sand to the glaciofluvial esker, which has a comparatively coarser texture and higher hydraulic conductivity.The concentration gradient at the front of the plume is steep in the esker (Fig. 11).The steepness can be explained by dilution, and by large differences in hydraulic conductivity between the esker and the surrounding material (Fig. 10).The distinct change from fine to coarse material alters the flow direction of the contaminated water to follow the border of the zone of coarse material.This affects the plume by narrowing the front as it moves towards the extraction wells.The limited lateral dispersal in the esker probably explains why the highly water soluble degradation product BAM can migrate in high concentrations over considerable distances.However, in comparison with an esker aquifer modeled by Artimo (2002) in Finland at similar scale, the longitudinal dispersivity values were 8 to 10 times smaller.Artimo did, however, not include the vadose zone in the modeling, and the aquifer was contaminated by the dense non-aqueous phase liquid tetrachloroethylene, which makes further comparisons difficult.
In a non-pumping scenario, the BAM-contaminated water plume still moves directly towards the extraction wells, in particular to the western well.In this well, the concentration increases when the pumps are stopped because the contaminated water becomes less diluted (Fig. 6).However, we can see that the opposite happens to the eastern well, where the concentration decreases when the pumps are stopped.According to the model the explanation is that the capture zone of the east well is partly located in finer and less permeable material (Fig. 10).Although we know that all samples since 2007 were taken in the western well, it is not certain which extraction well the older samples come from.If the results from the model are valid, it appears that samples older than 2003 were taken in the eastern well, but the remaining samples were taken in the western well.This assumption gives a correlation coefficient of 0.61 between observed and simulated concentrations, in comparison to 0.40 if all observations are assumed to come from the western well.
In order to speed up the cleaning process, the model was applied to optimize the location of remediation wells (Fig. 2).A 7-day pumping test was conducted in the newly established wells and the results validated the calibrated transmissivities and the predicted BAM concentrations.Transmissivities from the test were found to be 20 to 30 % higher than the calibrated model values (Table 5).The values from the pumping test may, however, have been overestimated because of intense rainfall during the pumping test (16 % of the annual net infiltration rate fell during that week).The predicted BAM concentration and the result from the pumping test agreed within the analysis error of 35 % for both wells.With the remediation wells, it was predicted, as a worst-case scenario, that the extraction wells supplying the municipal drinking-water utility would be clean in four years (concentration < 0.1 µg l −1 ).Without the remediation wells the same scenario predicted the extraction wells to be clean in nine years.
The observed and simulated concentrations were not perfectly matched.Reasons can be both that all significant transport processes were not fully understood and included in the model, and that the quality of the analyzed BAM samples in groundwater varies.In particular, the sampling quality may vary in those samples that were collected earlier than 2007.Since the model was run at regional scale combined with a deep vadose zone and coarse texture, it was assumed that preferential flow was insignificant.For example, it is not evident whether fingered flow can be a significant process at infiltration rates that normally occur at this site (Selker et al., 1992).However, any possible forms of preferential flow in the vadose zone were incorporated in the model through the dispersion term, using apparently higher macrodispersivity values than the true dispersion.In the groundwater model, the heterogeneous properties of glaciofluvial material were not included.According to the performed sensitivity analysis, the hydraulic conductivity turned out to have the highest effect on the model outcome.The variability of hydraulic conductivity may more or less influence the pathways of contaminant transport.If the variability is high, the assumption of homogeneous zones within the aquifer might have overestimated the predicted remediation efficiency (Maji and Sudicky, 2008;Rauber et al., 1998).
The distributed numerical approach applied in this study, can be compared to stationary, non-distributed model results described by Bergvall et al. (2007) for the same site.To model transport in the vadose zone, Bergvall et al. (2007) applied a simplified chromatographic flow model, implying that the infiltration rate and the total amount of water were constant.A retardation factor was used to account for sorption.To model transport in groundwater, Bergvall et al. (2007) used Darcy's law at five cross-sections along the flow direction between the nursery field and the two down-gradient extraction wells (Fig. 2).Dispersion was not directly modeled, but indirectly represented by the use of Monte-Carlo simulations to calculate contaminant transport within a range of values to give various outcomes.The degradation rate of dichlobenil was not represented in the model by Bergvall et al. (2007), and the results were highly dependent on sorption related coefficients.In groundwater, the transit times from the nursery to the extraction wells were found to be twice as long in the distributed model, compared to the simple nondistributed model.Mainly, this difference can be explained by differences in aquifer thickness and hydraulic conductivity.In particular, the thin layers of silt that partly cover the esker were not included in the non-distributed model.However, in the vadose zone the models gave similar median transit times for both water and BAM, with similar retardation factor of 1.7.This is probably a consequence of the thickness of the vadose zone, which results in a transit time of four years for water.Thus, the effect of the different resolution of the infiltration rate in the two models is small.An increased time resolution in the distributed model would be necessary to possibly obtain differences in median transit times.Despite the prediction ability of the simple mass-balance model in coarse-grained settings, the distributed numerical model should be used if the contaminant transport needs to be calculated more accurately, e.g. to predict concentrations at a specific point or to evaluate remediation strategies.
One limitation of our model is that breakthrough measurements of the contaminant at the extraction wells are not available, which may be crucial for predicting future transport behavior (Eggleston and Rojstaczer, 2000).However, the transport model is not only based on measured head and hydraulic conductivity values, but also on data from twelve years of contaminant measurements at several locations, including the source and the extraction wells.We believe that the results of this model show that it is possible to obtain satisfactory predictions of contaminants, despite the lack of breakthrough measurements, if supported by abundant field data.

Conclusions
This paper describes a coupled vadose zone-groundwater model for pesticide transport and its application to a pesticide-contaminated esker aquifer consisting of extensive glaciofluvial and littoral sediments.A sensitivity analysis revealed that hydraulic conductivity in the aquifer and infiltration rate account for almost half of the model uncertainty.At regional scale with the combination of a deep vadose zone and coarse texture, the observed concentrations could be described by the model without assuming preferential flow.However, a large macro-dispersivity was calibrated suggesting deviation from chromatographic transport in the vadose zone.The model approach is likely to be applicable at regional scale for similar aquifers since glaciofluvial and littoral sediments usually consist of well-sorted, unstructured sand and gravel.
The model was used to optimize the location of extraction wells for remediation.Running a worst-case scenario, the prediction showed that the establishment of two extraction wells would clean the aquifer in four years, compared to nine years without them.A pumping test of the two new extraction wells locally validated the calibrated transmissivities and the predicted BAM concentrations.
Further development of the model would require additional field measurements of macrodispersion to assess its importance in deep, sandy vadose zones.It is also important to characterize the variability of hydraulic conductivity and its effect on contaminant transport in eskers.

Fig. 1 .
Fig. 1.Modeling procedure and boundary conditions (BC) of water flow and solute transport.The vadose zone model was run only for the contaminated nursery area.

Fig. 2 .
Fig. 2. Location of the site within Sweden and map showing location of the nursery site in the modeled area, observation wells, extraction wells, deep soil samples, and salt injections.The contour lines of the hydraulic head are drawn every 5 m, and lakes, streams, erosion cuts, mires and exposed bedrock are represented.#7705 refers to an observation well; #1 and #2 to remediation wells.

Fig. 3 .
Fig. 3. Observed and simulated water contents in the soil profile that were sampled in 2003.The water content data obtained by the oven-drying method was used to calibrate the vadose zone model.The minimum, mean, and maximum values of water contents were obtained from measurements by the soil moisture probe at thirteen occasions (probe, min, mean, max), and was used for validation.The visualized mean of water content data obtained from the model (model, mean) refer to the corresponding thirteen days in the model.

Fig. 4 .
Fig. 4. Simulated and observed BAM concentrations at observation well #7705 located at the nursery field (see Fig. 2).The analysis error was ±35 % and is indicated by error bars.The moving average is the average of the closest four observations, before and after each observation.The observed peak value in 2001 was characterized as an outlier according to Dixon's Q-test.

a
Calculated by the Hazen method: K sat [m d −1 ] = 1000 × d 10 , where d 10 is the grainsize diameter [mm] at which 10 % by weight are finer.b Value refers to geometric mean.

Fig. 5 .
Fig. 5. Concentrations of BAM separated to central, western and peripheral plumes, in accordance with analyzed values (Table 2).The highest concentrations were found in the central plume, whereas the lowest concentrations were found in the peripheral plume.

Fig. 6 .
Fig. 6.Simulated and observed concentrations of BAM at the extraction wells.The analysis error was ±35 % and is indicated by error bars.Black bar at the bottom of the graph indicates that the wells were running; white bar indicates that they were closed.Note that before 2007, it is unclear whether the observed values were sampled in the western or eastern well.As from 2007, all observed values refer to the western well.The observed peak value in 2005 was characterized as an outlier, according to Dixon's Q-test.The European drinking water standard is 0.1 µg l −1 for pesticides.

Fig. 7 .
Fig. 7. Calculated infiltration rates (rainfall + snowmelt − evapotranspiration; left y-axis) and bottom water fluxes at 10.2 m depth (right y-axis).The daily water outflows were generated by the vadose zone model from 1 July 2003 to 30 June 2004.The indicated date of calibration refers to 1 November 2003.

Fig. 8 .
Fig. 8. Simulated and analyzed concentrations of dichlobenil at one identified hotspot (Fig. 2) in the vadose zone of the nursery field.The detection limit (d.l.) of the analyzed samples was about 10 µgl −1 .The date of calibration was 1 November 2003.The assumed application rate of dichlobenil was 1.7 g m −2 , applied at 24 times between 1976 and 1991.

Fig. 9 .
Fig. 9. Simulated mass fluxes (left y-axis) and concentrations (right y-axis) of BAM at the bottom of the 10.2 m deep vadose zone model.The daily outflows and concentrations of BAM were generated by the vadose zone model from 1 July 2003 to 30 June 2004.The assumed application rate of dichlobenil was 0.17 g m −2 , which was ten times lower than the rate in the calibrated model.The BAM concentrations from the vadose zone model were used as input to the groundwater model.

Fig. 10 .
Fig. 10.Boundary conditions of the groundwater flow model.The model consisted of two layers, which were divided into four texture groups with various hydraulic conductivities.The calibrated conductivities were as follows: 25 to 200 m d −1 for esker material, 2.0 to 20 m d −1 for sand, 0.0004 to 1.0 m d −1 for silt and 0.20 to 4.0 m d −1 for till.The modeled bottom elevation consisted mainly of the bedrock, but surfaces of silt or clay were found in some locations underlying more coarse-grained textures.

Table 3 .
Analyzed soil textures (refer to the soil profile sampled in 2003) and model parameters used in the vadose zone model: K sat is saturated hydraulic conductivity, χ is characteristic pore size, n is the tortuousity factor, T 1/2,dich and T 1/2 BAM are the half-lives of dichlobenil and BAM, respectively; and K d dich and K d BAM are the sorption coefficients of dichlobenil and BAM, respectively.main depth gravel a sand a silt a clay a K b vz

Table 4 .
Results of the sensitivity analysis for the concentration of BAM in the extraction wells.The sensitivities are expressed as percentages of the contribution to the total model uncertainty.The time step used in Eq. (5) was 180 days, and the time frame was 1 January 1976 to 31 December 2015.