Freely available online through the author-supported open access option.

Abstract

The spatial distribution of seepage at a flow-through lake in western Denmark was investigated at multiple scales with integrated use of a seepage meter, lake–groundwater gradients, stable isotope fractionation (δ18O), chlorofluorocarbon (CFC) apparent ages, land-based and off-shore geophysical surveys, and lake bed coring. Results were compared with a three-dimensional catchment-scale groundwater flow model using the MODFLOW and LAK3 codes for simulating lake–groundwater interaction. Seepage meter and model results of discharging groundwater to the lake compared well, if direct seepage measurements from near shore were combined with measurements from deeper parts of the lake. Discharge rates up to 9.1 × 10−7 m s−1 were found. Ground-penetrating radar used to map the lake bed sediments proved very effective in recognizing low- and high-permeability areas but also in understanding the complex recharge pattern of the lake and relating these to the geologic history of the lake. Recharge of the surrounding aquifer by lake water occurs offshore in a narrow zone, as measured from lake–groundwater gradients. A 33-m-deep δ18O profile at the recharge side shows a lake δ18O plume at depths that corroborates the interpretation of lake water recharging offshore and moving down gradient. Inclusion of lake bed heterogeneity in the model improved the comparison of simulated and observed discharge to the lake. The apparent age of the discharging groundwater to the lake was determined by CFCs, resulting in ages between 3 and 36 yr with an average of 16 yr. The simulated average groundwater age was 13.2 yr.

Seepage distribution at a flow-through lake was investigated using onshore and offshore geophysical, geochemical, and seepage measurements. Three-dimensional modeling results are in agreement with apparent groundwater ages and seepage rates. Offshore ground-penetrating radar combined with deep profiles of the natural abundance of oxygen-18 revealed an unexpected offshore recharge pattern from the lake.

Groundwater-dominated lakes are especially vulnerable to deterioration in lake water quality due to inputs from polluted groundwater. An understanding of the distribution and rate of seepage to and from lakes is therefore needed for environmental management or restoration of lake ecosystems (Hayashi and Rosenberry, 2002; Sophocleous, 2002; Gleeson et al., 2009). Several studies have demonstrated the influence of seepage on (i) lake water quality, e.g., the discharge of nutrient-rich groundwater (Loeb and Goldman, 1979; Brock et al., 1982; Belanger et al., 1985; Ito et al., 2007), groundwater rich in cations (Dean et al., 2003; Cullmann et al., 2006), dissolved inorganic and organic C (Striegl and Michmerhuizen, 1998; Staehr et al., 2010), or, in general, changes in lake alkalinity as a result of weathering processes in the watershed (Schafran and Driscoll, 1993); or (ii) biological communities, e.g., biodiversity and species distribution in seepage zones (Lodge et al., 1989; Hagerthey and Kerfoot, 1998; Rosenberry et al., 2000; Hayashi and Rosenberry, 2002; Sebestyen and Schneider, 2004).

The exchange of water and solutes between groundwater and lakes is complex and there is still a challenge in understanding the temporal and spatial variability across different scales (Käser et al., 2009). To address this challenge, a combination of monitoring and measurements at different scales for the whole lake and modeling is needed. Often, however, groundwater–lake interactions are only monitored or modeled for a small part of the lake and the spatial heterogeneity in seepage is not addressed at all. Kishel and Gerla (2002) found a high degree of seepage heterogeneity across just a 170-m2 area near the lakeshore. They concluded from a statistical analysis that even with 49 measurements of hydraulic conductivity, more data were required to cover the spatial heterogeneity in that relatively small area. Likewise, Guyonnet (1991) showed how small-scale sediment variations in relatively thin, either high- or low-permeability layers can have a significant effect on the seepage distribution. There is also contradicting evidence of the spatial distribution of seepage between lakes and groundwater in the near-shore environment both for flow from groundwater to lakes (herein referred to as discharge or inflow) and from lakes to groundwater (herein referred to as recharge or outflow).

Early theoretical and model studies predicted an exponential decrease in discharge with increasing distance from the shore (McBride and Pfannkuch, 1975; Winter, 1978); however, field studies have shown nonexponential relationships (Cherkauer and Nader, 1989; Shaw and Prepas, 1990; Schneider et al., 2005). On top of this, the regional hydrogeology also affects seepage. Model studies of whole-lake systems have been performed for hypothetical lakes, where system characteristics such as lake and aquifer geometry, slope gradients, and ratios of anisotropies in hydraulic conductivity (K) can easily be changed (e.g., Nield et al., 1994; Townley and Trefry, 2000; Smith and Townley, 2002; Turner and Townley, 2006). Before that series of studies, Winter (1983), Pfannkuch and Winter (1984), and Winter and Pfannkuch (1984) found that the factors needed to be known to understand lake–groundwater interaction are: (i) the water table configuration; (ii) the geometry of the groundwater system associated with the lake; (iii) the anisotropy of K; (iv) the distribution of high-permeability zones relative to the general geologic matrix; and (v) the lake depth and slope of the lake bed. Genereux and Bandopadhyay (2001) found from numerical studies that lake depths do not have a significant effect on the quantity or distribution of seepage. Furthermore they showed that adding low-conductivity lake sediments into a lake system has effects similar to increasing the anisotropy of the surrounding porous medium, both resulting in a evening out of the seepage rates throughout the lake bed. Cherkauer and Nader (1989) argued based on a field study that connectivity to an underlying bedrock aquifer controlled the seepage patterns. Schneider et al. (2005) found from comprehensive field studies that it was primarily the underlying geologic formations and regional precipitation patterns rather than lake bed sediments that controlled the observed seepage patterns. It is clear from these model and field studies that lake–groundwater exchange occurs at multiple scales.

Field studies tend to focus on seepage in the near-shore shallow zone, where accessibility is better than offshore. A wide range of field methods have been used, but data collection often represents a limited area of the lake. Conventional hydraulic gradient measurements based on either potentiomanometers or installations of piezometers are one of the most popular methods (Lee et al., 1980; Winter et al., 1988; Kishel and Gerla, 2002; Simpkins, 2006); seepage meters have been used in a number of studies (McBride and Pfannkuch, 1975; Lee, 1977; Cherkauer and Nader, 1989; Belanger and Kirkner,1994; Rosenberry et al. 2000; Kishel and Gerla, 2002; Schneider et al., 2005); tracers, such as the temperature in lake beds (Rosenberry et al., 2000; Kishel and Gerla, 2002), the temperature and conductivity along a lake bed (Lee, 1985), solutes (Lee et al., 1980; Frape and Patterson, 1981; Wollschläger et al., 2007), or plant community structure (Rosenberry et al., 2000) have been used to locate seepage zones.

The objective of this study was (i) to make a whole-lake evaluation and quantification of the seepage pattern at a flow-through lake by combining several different field methods with numerical modeling, (ii) to evaluate discharge and recharge fluxes at different scales, and (iii) to relate the observed seepage pattern to the geologic history of the lake. The study lake is Lake Hampen in the western part of Denmark. It has received some attention from ecological studies in the past (Christiansen et al., 1985; Christensen and Sorensen, 1986), but the hydrology and hydrogeology of the lake is largely unknown. Moeslund (2000) estimated the water budget of Lake Hampen and concluded that it was possible that the lake received a large inflow of groundwater but did not have any field data to support this statement. A suite of methods were used at different scales in the current study: direct (seepage meter) and indirect methods (tracers, δ18O and CFC ages, vertical lake–aquifer hydraulic gradients, catchment well–lake gradients, lake bed geologic mapping, and geophysical methods). Field measurements were compared with a catchment-scale numerical model, where lake–groundwater interaction was simulated using the MODFLOW LAK3 lake package (Merritt and Konikow, 2000). Special emphasis has been on also addressing the recharge patterns at the lake. The literature on seepage and flow-through lakes is abundant, but despite this, recharge from lakes has seldom been described in field studies with only a few exceptions (Rosenberry, 2000, 2005).

Site Description

Lake Hampen is situated in the middle of the Jutland Peninsula in the western part of Denmark and forms the upper portion of the Skjern River catchment (Fig. 1 ). The lake has a surface area of 760,000 m2 and a volume of 3200 ML. The lake consists of a western shallow basin and an eastern deep basin (Fig. 2a and 2b ) with a maximum depth of 13.2 m. The catchment is largely covered by forests and grass fields without any urban areas.

The regional geology of the area is defined by the main advance of the last glaciation in Denmark or the Late Weichselian glacial maximum (Houmark-Nielsen, 1989), which occurred 18,000 yr ago according to Larsen and Kronborg (1994). In the middle and southern parts of Jutland, the main advance is bisected by a north–south-trending line, with glacial outwash situated west of the line and clayey tills east of the line. Lake Hampen is located east of the border (Fig. 1) and is a part of the so-called Lake Highland in the central part of Jutland (Larsen and Kronborg, 1994). Lake Hampen was probably formed as a dead ice hole when the ice cover retreated toward the east, leaving a 15- to 35-m layer of coarse meltwater sand and gravel in the area close to the lake. This can be observed at two sites, Multilevel Well (MLW) 4 and MLW3, 10 and 50 m east and west of the lake, respectively (Fig. 2a and 2b) and from geophysical surveys (Fig. 3 ). Tertiary marine silty clays are present in MLW4 and MLW3 below the coarse glacial sediments. Between the upper glacial sand and the lower quartz sand there is a hydraulic head difference of 4.5 m separating the two layers into an upper and a lower aquifer.

The hydraulic conductivity of the upper aquifer near the lakeshore at the eastern side of the lake (4.3 × 10−3 m s−1 at MLW4) is somewhat higher than at the western side (2.6 × 10−3 m s−1 at MLW3). Farther to the west of the lake, however, the hydraulic conductivity decreases to 1.8 × 10−4 and 5.2 × 10−4 m s−1 at Boreholes B7 and B8, respectively (Fig. 2a). This pattern could be explained by the depositional environment after the ice cover retreated from the lake area, with the glacial front located just east of the lake (Larsen and Kronborg, 1994). The coarsest sediments would be deposited near the glacial front and sediments would become finer at greater distances from the glacier.

The lake sediments consist of fine-grained organic sediments (predominantly gyttja) in the deeper parts of the lake (Moeslund, 2000). In the western shallow part of the lake, a peat layer has been observed and interpreted to originate from a period with a lower water stage (Moeslund, 2000). This peat layer possibly exists throughout the shallow western basin. Because of recreational use, sand has been dumped along the western shoreline, covering the lake bed from the shoreline and extending presumably 15 to 50 m offshore.

Precipitation measured at a climate station 1.5 km south of Lake Hampen averaged 888 mm yr−1 for the period 1970 to 2007, with a calculated standard deviation for the period of 143 mm (Danish Meteorological Institute Station 23030). During the study period, February 2007 to December 2009, the lake was only ice covered from 5 to 19 January 2008 and even then, smaller ice-free areas could be found. The calculated potential evaporation from the lake shows an annual evaporation from the free lake surface of 407 mm in 2008. This can be compared with data from Moeslund (2000), who estimated evaporation from the lake in 1998 to be 425 mm.

Lake Hampen has two minor inlets in the eastern part of the lake and a minor outlet in the western part (Fig. 2a). The two inlets supply water from two small lakes or wetlands, Kragsoe and Troldsoe. During the study period, Troldsoe was dammed by the forest authorities and the inlet was dry. The inlet from Kragsoe was observed monthly to detect any increase in flow but it did not exceed 0.5 L s−1. The outlet flow was measured continuously using a V-notch weir 150 m downstream from the lake, causing a minimal backwater effect, and just upstream of the weir. The backwater effect was kept between 5 and 15 cm at the weir (upstream side) by changing the angle. During 2008, the outlet flow averaged 4.8 L s−1 with maximum flow of 15.6 L s−1 in early April. The outlet dried out 1 wk in August. Moeslund (2000) reported that in 1998 the outlet was dry throughout the year, but 1.2 and 0.1 L s−1 entered the lake through the inlets from Kragsoe and Troldsoe. Hydrologic variables are summarized in Table 1 .

Experimental Methods

Climate Station

A climate station was placed ∼35 m offshore at the lake (Fig. 2b). The station measured temperature, precipitation, and incoming and outgoing short- and long-wave radiation. Precipitation was measured with a tipping bucket rain gauge with a built-in datalogger (Eijkelkamp Agrisearch Equipment, Giesbeek, the Netherlands). Precipitation was corrected with monthly factors recommended for Danish conditions by Allerup and Madsen (1979). The radiometer was a Kipp & Zonen CNR 1 (Kipp & Zonen, Delft, the Netherlands) with a spectral range of 300 to 2800 nm. The resistance-based Penman–Monteith equation (Shuttleworth, 1992) was used to calculate daily evapotranspiration rates for conifer forest and grassland. Evaporation from the lake surface was calculated according to recommendations made by Shuttleworth (1992) for calculation of the potential evaporation from free water surfaces. Groundwater recharge was estimated on the basis of four different scenarios: recharge in forested or grass areas through a 2- or 5-m unsaturated zone. The HYDRUS-1D version 4.06 code (Šimůnek et al., 2008) was used to simulate recharge to the aquifer. Evapotranspiration estimated by the Penman–Monteith equation was used as input to the HYDRUS-1D model. The Feddes root uptake model was used with default values from the HYDRUS-1D database for the actual vegetation (Šimůnek et al., 2008).

Potentiomanometer and Piezometer Transects

A hydraulic potentiomanometer, similar in principle to the one described by Winter et al. (1988), was used to establish the hydraulic gradients near the lake perimeter during October 2006 to March 2007. During that preliminary survey of lake–aquifer gradients, the small screen of the potentiomanometer, 3.5 cm in length and with an outer diameter of 1.7 cm, often was clogged by fine-grained particles. This was experienced in the shallow western part of the lake even though the top layer was sand. During 2 to 3 Sept. 2009, the gradient measurements were repeated, but this time with steel piezometers with 9-cm screens and an outer diameter of 2.8 cm. The steel piezometers were forced 50 cm below the lake bed at a water depth of 50 cm and the difference between the hydraulic head at the piezometer screen and the lake stage was measured after the hydraulic head in the piezometer stabilized. During this survey, it was observed by using slug tests (falling head) that the near-shore permeability in the lake bed was considerably lower in the shallow basin than in the deep basin (head recovery times were relatively long, 30–60 min, in the sediments below the shallow basin compared with recovery times of a few minutes in sediments below the deep basin).

Five piezometer stations (P1– P5), with one to three piezometers at each station, were installed around the lake (Fig. 2b). The polyvinyl chloride (PVC) piezometers have an outer diameter of 2.7 cm and a screen length of 10 cm. The piezometers were installed by a motorized auger, drilling a 7.5-cm-wide borehole. The remaining borehole space was filled with medium size sand in the filter section, with a bentonite seal on top. The auger material was backfilled in the rest of the borehole. The piezometers were placed 2 to 50 m from the lake shoreline.

Eight piezometers (B1–B8) were installed in the catchment surrounding the lake 1.5 m below the groundwater table with a GeoProbe (GeoProbe Systems, Salina, KS). These piezometers were made of PVC with an outer diameter of 3.2 mm and a screen length of 1 m. The purpose of B1 to B8 was to gather additional data on hydraulic heads, supporting the existing wells placed in the upper aquifer (B10–B21, Fig. 2a).

Based on the observed hydraulic gradients in the near-shore environment, two transects were established in a discharge zone (I2 and I3) and one transect in a recharge zone (O2) to investigate seepage at the local scale. The transects were also used to investigate the temporal hydraulic gradient patterns between the upper aquifer and the lake. The three transects included 12 to 15 piezometers. Of these, three to four piezometers were installed in the lake and eight to 12 on land with a distance to the shoreline of up to 70 m. The piezometers installed in the lake were made of 2.8-cm-diameter galvanized steel pipes with a 9-cm-long steel screen and installed with a pneumatic hammer. The screens of the transect piezometers were placed 0.5 to 1.0 m below the groundwater table or lake bed. The type and installation of the piezometers on land were similar to the P1 to P5 wells. Hydraulic heads were measured monthly in 2008. At the same time, hydraulic heads were measured in the 10 piezometers at sites P1 through P5, at boreholes B1 to B8, and at the nine existing wells (B9–B21) in the catchment (Fig. 2a and 2b).

In July 2008, steel piezometers were installed offshore in the western shallow part of the lake and offshore near Transects I2 and I3 in deeper parts of the lake (Fig. 2b). The piezometers were removed again after measuring the vertical hydraulic gradient, taking water samples, and coring the upper approximately 1 m of lake bed with a half-core Russian sampler.

Water Samples

Analysis of δ18O was done according to the convention described by Appelo and Postma (2005): 
\[\mathrm{{\delta}}^{18}\mathrm{O}_{\mathrm{sample}}{=}\frac{(18_{\mathrm{O}}/16_{\mathrm{O}})_{\mathrm{sample}}{-}(18_{\mathrm{O}}/16_{\mathrm{O}})_{\mathrm{VSMOW}}}{(18_{\mathrm{O}}/16_{\mathrm{O}})_{\mathrm{VSMOW}}}1000\]
[1]
with Vienna Standard Mean Ocean Water (VSMOW) having a 18O/16O ratio of 2.005 × 10−3 (Appelo and Postma, 2005). The δ18O samples were analyzed with a dual-inlet Isoprime mass spectrometer (Isoprime Ltd, Cheatle Hulme, UK) with a Gilson Multiflow unit and results given as per mil (‰). Samples were analyzed five times and a standard deviation was calculated for each sample. Water samples used for the analysis of δ18O were taken with a vacuum pump and kept in 20-mL polyethylene bottles at 5°C until analysis. The water volume inside the piezometers was replaced a minimum of five times before sampling. Water sampling was done at MLW1 and MLW2 (Fig. 2) in July 2008 using the described steel piezometers driven down with a pneumatic hammer. The deepest samples were taken at 9 and 11 m below the surface. The steel pipes had lengths of 2 m and after every second sample a new 2-m section of pipe had to be fitted through a threaded pipe-connector fitting. The threaded pipe-connector fitting had a 6-mm-wider outer diameter than the steel pipes. The procedure was also attempted at a third site near Transect O2, but after taking the second sample at 2 m, the screen clogged and sampling below 2 m was not possible. This is probably because of an impermeable horizon near the 2-m depth. Water sampling was done every 1 m down to 18 and 33 m at MLW4 and MLW3, respectively, using a GeoProbe with a 30.5-cm-long screen and an outer diameter of 2.5 cm. A peristaltic pump was used to sample the water from the 12-mm polyethylene tube fitted to the screen. Hydraulic tests and sediment cores were also taken at MLW3 and MLW4.

Water samples used for chlorofluorocarbon (CFC) analysis were taken at MLW4 in flame-sealed ampoules, using the sampling procedure described by Busenberg and Plummer (1992). A second set of samples was collected in screw-cap bottles. A 100-mL bottle placed in a 1-L container was filled from the bottom, allowing the steady flow of water to fill the container before closing the bottle under water. The samples were taken to the laboratory and analyzed according to the procedure described by Busenberg and Plummer (1992). Field measurements of dissolved O2 content and pH were performed along with groundwater sampling.

Seepage Meters

Seepage meters were constructed as described by Lee (1977) to measure seepage through the lake bed. The half-barrel steel seepage meter had an inner diameter of 57.5 cm and a height of 25 cm. A steel pipe with an inner diameter of 2.7 cm was fitted on the top of the seepage meter, with a Storz coupling at the opposite end of the pipe. The Storz coupling is a type of hose coupling composed of interlocking hooks and flanges. The Storz fitting is used to quickly connect two pipes and is made of aluminum with a rubber gasket between the two opposite parts of the fitting. The seepage meter was installed with the top of the half barrel 5 cm elevated from the lake bed to ensure free flow toward or from the pipe inside the meter. A plastic bag was attached to a 3-cm-long piece of steel pipe with a rubber band. The diameter of the steel pipe was 2.7 cm. On the other end of the pipe, a two-way valve with a Storz fitting was fitted. After prefilling with 1 L of water, the plastic bag arrangement was attached to the equilibrated seepage meter and the valve was opened. A 20- by 40-cm transparent Plexiglas shield was gently placed over the plastic bag and attached to the meter to hold the bag in place and minimize velocity-head effects associated with waves or currents flowing past the plastic bag (Rosenberry, 2008). The measurement time was from 30 min to 10 h depending on the seepage rate.

Previous uses of seepage meters (Erickson, 1981; Cherkauer and McBride, 1988; Belanger and Montgomery, 1992) have indicated that the measured seepage meter rates should be corrected using a correction factor. Rosenberry and Menheer (2006) reviewed the use of correction factors and found reported values between 1.05 and 1.74. Belanger and Montgomery (1992) reported a ratio of 0.77 between the seepage meter value and the actual seepage at fluxes below 6 × 10−6 m s−1. This gives a correction factor of 1.23, which was applied to the measurements in this study because the seepage meters and procedure used by Belanger and Montgomery (1992) were similar to the ones used in this study (except that the inner diameter of the tubing connecting the plastic bag to the meter was smaller in Belanger and Montgomery [1992] at 5.6 mm). All of the measured seepage meter rates were below 6 × 10−6 m s−1 at Lake Hampen. Seepage meter measurements were performed within 31 m of the shoreline and at water depths below 1.2 m, except at I3 where seepage meters were installed at water depths of 4 and 6 m.

Geophysical Investigation

The geology near Lake Hampen was examined by 11 multielectrode profiles (MEPs; Dahlin, 1996) in the vicinity of the lake (Fig. 2a). The survey used SYSCAL PRO equipment (Terraplus USA, Littleton, CO) with 84 electrodes at 5-m spacing. A Wenner measurement scheme, with 756 measurements, was used, giving a penetration depth of almost 70 m in the central section of the profile and a total length of 420 m. The Res2Dinv software program (Loke, 2008) was used for inversion of the measured data. The profiles were interpreted using the typical resistivity values given by Ahrentzen et al. (1987) for Danish soils and geologic units.

Natural γlogging of B9 and B10 (Fig. 2a) was done with a Robertson Geologging probe (Probe no. 25 061 000, Robertson Geologging, Deganwy, UK) using a 50- by 25-mm scintillation crystal as a natural γ detector. ViewLog 3.0 software (Earthfx Inc., Toronto, ON, Canada) was used for plotting the measurement data.

Deposition of fine-grained organic deposits can greatly affect the distribution of seepage in lakes (Winter et al., 1998) and therefore a ground-penetrating radar (GPR) survey was performed during August 2008 to detect the presence and thickness of organic sediments in the lake. The technique of using GPR to survey in-lake geology was described in studies by Caldwell and Fitzgerald (1995) and Buynevich and Fitzgerald (2003). The GPR lines were collected using a PulseEKKO Pro System (Sensors & Software, Mississauga, ON, Canada) equipped with two 100-MHz antennas. The water depth was calculated assuming an electromagnetic wave velocity of 0.033 m ns−1. The antennas were placed in a separate inflatable boat with a 1.5-mm-thick rubber bottom to minimize signal attenuation. A total of 32 lines was collected. The GPR measurements were correlated to the sediment cores taken with the Russian sampler in July 2008.

Numerical Methods

A numerical, three-dimensional, steady-state model using the MODFLOW code (McDonald and Harbaugh, 1988) was set up for Lake Hampen and the catchment surrounding the lake. The LAK3 lake package (Merritt and Konikow, 2000) was used to simulate the interaction between Lake Hampen and the surrounding aquifer. The LAK3 package was evaluated by Hunt et al. (2003) and found appropriate to simulate lake–groundwater interaction. The northern, eastern, and southern model boundary was defined by the groundwater divide, shown as the catchment boundary in Fig. 1. The western model boundary was defined with a specified head based on data from boreholes along this perimeter (Fig. 4a ). Recharge zones are shown in Fig. 4b, where Zones 1 and 2 are forest and grass cover, respectively, with a 0- to 5-m unsaturated zone, and Zones 3 and 4 are forest and grass cover, respectively, both with an unsaturated zone >5 m thick. The LAK3 package requires inputs of precipitation, calculated potential evaporation from the lake surface, net flow to or from the lake via surface water, and a parameter defining the lake bed permeability (hydraulic conductivity divided by lake bed thickness). The geology was represented in the model as four hydrogeologic units: (i) glacial sand, the upper aquifer 15 to 35 m thick; (ii) glacial clay, the top layer in the most northern part of the catchment and in some areas between the glacial sand and tertiary clay; (iii) clay, from 10 to 35 m in thickness below the glacial sand or clay; and (iv) sand that forms the bottom layer in the model. Horizontal discretization was 100 by 100 m, but with a finer 50- by 50-m grid near the lake (Fig. 4a). A 10- by 10-m digital elevation model was used to interpolate the model top boundary. In the same way, a 10- by 10-m digitized bathymetrical map of Lake Hampen was used to interpolate the lake bed elevation to the 50- by 50-m model grid. The model had eight numerical layers.

The model represents steady-state conditions in the year 2008. This means that every input to the model was calculated as average values in 2008. The model was calibrated against 32 head measurements, all located in the upper aquifer and the lake surface (Fig. 4a) (11 measurements were made at each of the 32 sites during 2008). The LAK3 package calculates a separate water budget for the lake and a lake stage. Standard error indicators such as mean error (ME), mean absolute error (MAE), and root mean square error (RMSE) using computed and observed heads were used in the calibration assessment. Emphasis in the calibration process was on the model's ability to match the observed mean lake stage of 79.06 m above sea level.

Vertical anisotropy in hydraulic conductivity was kept constant during the calibration with a value of 3. A vertical anisotropy of 2 was used by Hinsby et al. (2007) for a similar sandy aquifer in Jutland. Horizontal isotropy within each unit was assumed.

Two calibrations were performed, one where the lake bed had a homogeneous leakage parameter (Merritt and Konikow, 2000) for the entire lake and one where two leakage parameters were defined. In the second calibration, areas covered by the fine-grained organic deposits or the reported peat layer in the shallow part of the lake were assigned a very low leakage parameter equivalent to a hydraulic conductivity of 3 × 10−9 m s−1. The leakage parameter for the remaining more permeable areas near the shoreline was calibrated. For convenience, the two calibrations are referred to as the homogeneous and heterogeneous models, respectively. After calibration, the model was used to compare seepage through the lake bottom with the observed seepage meter measurements. The MODPATH code (Pollock, 1994) was use to estimate the age of the groundwater discharging to the lake for the heterogeneous model. One hundred particles were tracked backward from cells in Layers 1 to 5 at the location of MLW4. A porosity of 0.3 was assumed.

Results

Hydrologic Parameters

Precipitation measured at Lake Hampen was 901 mm in 2008, which is close to the mean value of 888 mm yr−1 observed between 1970 and 2007 (Danish Meteorological Institute Station 23030). The calculated recharge for the four different zones varied between 412 and 622 mm (Table 1). The difference of 210 mm yr−1 between coniferous forest with an unsaturated zone of 5 m and grass cover with an unsaturated zone of >5 m was expected because of the possible high interception loss in coniferous forests. Zones 1 and 2 show a difference of 161 mm between forest and grass cover with an unsaturated zone of 0 to 5 m. Zones 3 and 4 show a difference of 140 mm with an unsaturated zone > 5 m. These results are comparable to a study made by Stewart (1977), where the net additional evaporation due to interception was 145 mm for a pine forest. No significant surface inflows were observed during 2008 from the small channels from Kragsoe and Troldsoe (Fig. 2a). The measured outflow during 2008 at the outlet V-notch weir was 1.52 × 105 m3 yr−1, equivalent to 200 mm yr−1 normalized with the lake area.

Hydraulic Gradients and Stable Isotope Fractionation

Transience in hydraulic heads measured in wells near the lake and in the catchment show similar patterns, with the highest values in April and decreasing toward September and October, when significant recharge to the aquifer starts again (Fig. 5 ). The lake water stage also followed the yearly pattern observed in the aquifer (Fig. 5). The yearly variation in hydraulic head in the upper aquifer was greater the farther the distance from the lake. The amplitude in hydraulic heads in B11 (for location see Fig. 2a) was 1.07 m, whereas the amplitude at B6 was only 0.31 m. The water stage in Lake Hampen changed 0.29 m during the study. The hydrograph of B2 is somewhat different from the others. It seems that the yearly rising and falling of piezometric heads is several months delayed. One possible explanation is that the unsaturated zone of ∼10 m causes a delayed piezometric response in comparison with the other wells. Second, Borehole B2 and B3 are close to the border of the interpolated groundwater divide, which may be affected by errors. The precise determination of the groundwater divide was difficult with only a few wells in this area.

Vertical hydraulic gradients in the near-shore zone at Lake Hampen are shown in Fig. 6 and Table 2 . A hinge line, as defined by Winter (1999), dividing the lake into a discharge and a recharge area, is difficult to establish. Generally, positive gradients were observed in the northeastern part and negative gradients in the southwestern part of the lake. High gradients up to 0.12 were observed near Transects I2 and I3 and in the southern part of the lake, but the majority of positive gradients were between 0.01 and 0.03. The majority of negative vertical gradients had values down to −0.045, with a single observation down to −0.208 (Fig. 6).

Vertical hydraulic gradients observed during the offshore campaign, July 2008, in the shallow western part of the lake, showed high negative gradients of −0.421 to −0.331 at P6 to P9 (Fig. 2b). The lake bed at P6 to P9 had an upper layer of fine-grained gyttja and highly decomposed peat that ranged in thickness from 65 to 117 cm. The screens at P6 to P9 were installed below this horizon. At P10, the hydraulic gradient was only −0.086 and the lake bed consisted of gravel and coarse-grained sand. Offshore measurements at I2 and I3 showed positive gradients of 0.009 to 0.033 and resembled the ones in the near-shore zone (Table 2).

The δ18O in Lake Hampen varied between −3.25 and −4.47 and differed from groundwater samples taken in catchment wells, which had an average value of −8.25 (Table 2). The calculated standard deviation among the seven aquifer samples was 0.15. Values from Transects I2 and I3 clearly indicated groundwater composition (Table 2; Fig. 7 ). Transects I2 and O2 (Fig. 7) show the difference between discharge and recharge zones, with values close to groundwater at I2 and values ranging from the lake water signal toward a mixed signal farther inland. At O2, values greater than lake water values were observed, indicating a further fractionation within the lake bed or lake water conditions, with δ18O values less than measured.

Figure 8 shows the vertical profiles of δ18O in the discharge (MLW4) and recharge areas (MLW1, MLW2, and MLW3). Multilevel Well 4 located near I3 showed values between −8.83 and −7.22, with an average of −8.14 for the 16 samples taken from the 2- to 17-m depth. The data show a slightly decreasing trend with depth. The average value for the vertical profile is close to the −8.25 observed for the catchment samples. Transect I3, with values of −10.65 to −6.94, is located close to MLW4. The average δ18O value for I3 is −7.72, which is comparable to the upper part of MLW4.

The three vertical profiles in the recharge zone, MLW1, MLW2, and MLW3 (Fig. 8), generally showed a similar pattern, with a close-to-groundwater signal of −7 down to the 6- to 8-m depth and then a change toward a mixed groundwater–lake water signal. Between 12 to 15 and 20 to 26 m, MLW3 showed values above −4.02, indicating the presence of lake water. The deepest observation point in the profile was at 33.5 m with a value of −7.29 close to that observed in the groundwater elsewhere. The MLW1 and MLW3 wells are only separated by ∼20 m and the two profiles show a very similar pattern for the depths covered by MLW1, although samples at MLW1 and MLW4 were taken in July 2008 and those at MLW3 in July 2009. Values down to 7 m in MLW2 (−7.08 to −6.19) are similar to the ones at P2.1 to P2.3 (−7.22 to −5.47) and indicate a groundwater signal with some evaporative fractionation or mixing with lake water.

Geophysical Measurements

The objective of the land-based geophysical investigation was primarily to support the conceptual geologic model of the Lake Hampen catchment used in the modeling work. Figure 2a shows the 11 MEP lines and Fig. 3 illustrates the output from a selected MEP investigation. The profile is interpreted with an upper unsaturated zone of glacial meltwater sand ranging from 0 to 10 m below the surface, with resistivity values >1000 Ω m. Below this, the glacial sand became saturated and extended down to 25 to 35 m, with resistivity values of 120 to 1000 Ω m. The lower part of the profile was interpreted to be clay with a high content of silt in the upper horizons of the layer because normal resistivity values for Danish clay would be below 70 to 80 Ω m (Ahrentzen et al., 1987). At B9 (Fig. 2a), the clay was observed to be 41 m thick and a γ log showed that the silty clay observed at MLW4 was becoming progressively more fine-grained with depth. The MEP lines indicate that the silty clay layer is persistent throughout the catchment. Below the clay, at least 14 m of quartz sand was observed at B9. Glacial clay (till) was observed between the glacial sand and the clay in some areas of the catchment but also as the surface layer in the most northern part of the catchment.

The 32 GPR lines resulted in a map of the lake bed, where thicknesses of the organic layer in most parts of the lake were estimated. Figure 9a shows GPR Line 31 starting in the shallow western part of the lake and moving eastward (for location, see Fig. 2b). The red circles show where the organic layer seems to terminate or overlaps the underlying sandy aquifer or beach sediments. The fine-grained organic layer was interpreted to be between I and II, in the shallow part of the lake, and between III and IV in the deep part and is marked with green arrows. Furthermore, a thin layer seems to extend to the western shore from the threshold. Following the uppermost reflector (white line, illustrated by the blue arrow) from position 0 m toward the east, the reflector overlaps at II (western side of the threshold). This indicates that there also exists an upper organic layer (peat) that extends onto the beach and does not terminate at Position I. From Fig. 9a, it can be seen that deposition of fine-grained organic matter started at shallower water depths at Position III than at IV. This is because the general wind direction in Denmark is westerly. The GPR Line 12 is located in an area of the lake where the lake bed slopes steeply (Fig. 2b). Because of the water depth and the thickness of the organic layer, the profile does not penetrate the organic layer and reflection lines do not extend beyond about 5 m beneath the lake bed (Fig. 9b). Areas with an organic layer thickness <5 m are shown in Fig. 9c.

Seepage Meter, Chlorofluorocarbons, and Modeling Results

The two calibrations of the groundwater models resulted in similar error values and acceptable precision in the calculation of the lake water stage (Table 3 ). Computed heads were nearly identical for the two calibrations and close to the observed values (Fig. 10 ). Hydraulic conductivities found by calibration were within normal ranges for Danish geologic units defined by Henriksen et al. (2003).

The homogeneous and the heterogeneous calibrations resulted in leakage parameters of 3.0 × 10−6 and 6.1 × 10−6 s−1, respectively. The difference is a result of the reduced lake bed area in the heterogeneous model, where exchange between the lake and the aquifer occurred in a smaller area.

Figure 11 shows the modeled shoreline seepage for the two models, together with the results from seepage meter measurements at the shoreline. Peak discharge and recharge for the homogeneous model were generally lower than for the heterogeneous model. Differences between the two models are primarily seen for the recharge pattern at the threshold between the two basins and in the shallow part of the lake. The homogeneous model had peak recharge in the most western part of the shallow basin, whereas seepage was shifted toward the threshold in the heterogeneous model. The modeled seepage from the heterogeneous model (in the area where it flattens to zero in Fig. 11) should not be compared with measurements because it is a result of the imposed impermeable leakage parameter in the heterogeneous model. In terms of discharge, the striking difference between the two models is that the heterogeneous model simulated more extreme high and low discharges. The discharges measured at I2 and I3 were simulated slightly better by the more extreme high and low pattern of the heterogeneous model.

Concentrations of CFCs from MLW4 near the eastern shoreline (Fig. 2a) is shown in Table 4 together with the apparent CFC ages calculated by the method of Busenberg and Plummer (1992). The CFC partitioning between soil, air, and groundwater was assumed to have happened at 8°C, the annual mean temperature in Denmark. As a result of the Montreal protocol concerning the phasing out of the production of CFCs, the atmospheric CFC concentrations are now slowly decreasing, after having reached their maximum in 1994 (CFC-11 and CFC-113) and 2002 (CFC-12). Therefore, theoretically the minimum age (maximum concentration) to be calculated using the CFC method in 2009 is 15 and 7 yr, respectively. In practice, however, when taking into account the variation in analytical results of repeated samples (1–5%), the minimum age to be determined with any reasonably confidence is approximately 10 yr. Furthermore, two different CFC ages may be assigned to certain CFC-11 concentrations across a narrow range due to the decrease in atmospheric CFC-11 in recent years (Table 4).

Contrary to what might have been expected, dissolved O2 gradually increased with depth from 0.4 mg L−1 at 2 m, just 1.5 m below the water table, to 8 mg L−1 at 9 m (Table 4). Nitrite (NO2) and N2O were present in the upper 5 m, indicating that denitrification occurs at shallow depths but not deeper. Nitrous oxide was detected during the CFC analyses but not quantified because it appeared as peaks in the chromatographic curve, as described by Busenberg and Plummer (1992). According to Hinsby et al. (2007), partial degradation of CFCs, particularly CFC-11, may occur during denitrification, leading to higher CFC ages. This may be one reason for the higher CFC age based on CFC-11 relative to the other two CFCs in the upper 5 m. To be able to compare groundwater residence time with that predicted by modeling a single CFC, age has been assigned to each depth of sampling (Fig. 12 ). For the upper 11 m, the CFC age was primarily based on CFC-12 results, while for deeper levels both CFC-11 and CFC-12 were used.

Besides dissolved O2, other solutes that showed a significant change with depth include K, Ca, NO3, and HCO3, the first three of which are shown in Table 4. The changes observed in Ca and NO3 may fairly easily be explained as being due to the dissolution of CaCO3 in different soils and the application of fertilizers. The extremely high K concentration, however, up to 110 mg L−1 in the upper 6 m, must come from a local source, most likely the old farm 30 to 60 m upstream of the MLW4 location.

Results from the particle simulations are shown together with CFC results in Fig. 12.

A frequency distribution for Layers 1 to 5 is shown in Fig. 12b. Layer 1 was simulated to receive relatively young water of 0.5 yr. Layers 2 to 5 show a two-peak distribution. Particles arriving in Layer 2 have ages between 1 and 16 yr and peak frequencies of 18 and 6% with ages of 6 and 15 yr, respectively. Layer 3 particles have ages of 6 to 20 yr and peak frequencies of 16 and 7% with ages of 10 and 18 yr, respectively, but 57% of the particles in this layer have ages between 7 and 10 yr. Layer 4 also shows a two-peak frequency distribution of 14% with ages of 14 and 15 yr and 5% with ages of 21 yr. Layer 5 shows a relatively flat distribution, with ages ranging from 13 to 53 yr. The mean age for Layer 5 is 24.9 yr.

Discussion

The observed patterns in hydraulic head in the catchment and lake water stage throughout 2008 indicate a good hydraulic connection between Lake Hampen and the surrounding upper aquifer. Because the lake water stage and hydraulic heads in the surrounding aquifer vary synchronously, the horizontal hydraulic gradient between the lake and the aquifer appears to be close to constant. This does not imply that seepage is constant at a given portion of the lake bed but indicates that the assumption of modeling the system as a steady state is reasonable (Fig. 5). The relatively small temporal variations in the gradient between the aquifer and the lake water stage suggest that great temporal changes in seepage patterns are not the case at Lake Hampen.

Water sampling in a system like a groundwater–surface water system is always associated with some uncertainty because the water sample represents only a momentarily image of the system. The water sampling was not done simultaneously at all the different scales involved, which therefore adds uncertainty in the interpretation. It is believed, however, that the stable seepage pattern is representative for yearly system variations, and therefore water samples collected at various times can still be used for analysis of the system.

The general trend of δ18O values in the catchment and at MLW4 clearly shows that groundwater discharge to the lake has values that differ from the lake signal (Table 2; Fig. 8).

Samples taken at MLW3, west of the lake, indicate a very deep lake water plume with an unmixed lake signal between 12 and 15 and between 20 and 26 m. This is supported by MLW1 and MLW2 for the upper 9 to 10 m, even though they were sampled a year before MLW3.

The geophysical investigation and the geologic mapping of the lake bed in the shallow western part of the lake indicate a peat layer, which terminates at the threshold between the deep and shallow basins. The peat layer was probably formed during a period of the lake's history with a lower water stage (Moeslund, 2000). It is very possible that the wetland generating this peat layer originally extended to the recent perimeter of the shallow part of the lake and that the relatively low permeability observed is a result of the decomposed peat that accumulated during the period of lower water stage. The consequence will be that the zone where lake water recharges to the aquifer is shifted offshore to the margin of the peat deposit. The vertical profile of δ18O at MLW3 seems to confirm this recharge pattern because the unmixed lake signal appears at the deeper levels in the profile (Fig. 8). Another evidence for this offshore recharge pattern is the vertical gradients observed between the lake and the sand horizon just below the peat layer seen at P6 to P10. Piezometers P6 to P9 (Fig. 2) all showed a very high negative gradient between −0.331 and −0.421, but just beyond the eastern margin of the peat deposit, at P10, the gradient suddenly decreased to −0.086.

If no significant exchange of lake water to the aquifer occurs west of the threshold in the shallow basin, the observed values in the upper part of the MLW1, MLW2, and MLW3 profiles would be unfractionated recharged precipitation (originating from recharging precipitation between the lake shore and the boreholes). The upper δ18O values of MLW1, MLW2, and MLW3 and values from P2.1, P2.2, and P2.3 are slightly higher than the groundwater observed in the catchment. The above-groundwater levels could be caused by flooding with the lake extracting and directly diluting the precipitated recharge with lake water. A seepage face at the new shoreline would then be possible in the new areas not underlain by the decomposed peat layer, enabling high recharge rates. This could lead to the mixed signal observed in the upper parts of the profile. The scenario is similar to the one proposed by Rosenberry (2000) where the normal lake shoreline is exceeded during high water stage with a resulting rapid seepage face because the new lake area is not underlain by low-permeability lake sediments. Another explanation could be that diffusion mixes the lake and rainfall signals to some degree. A third, and maybe the most obvious, explanation is that the low seepage rates occurring near the western shoreline slightly dilute the otherwise undiluted infiltrated precipitation. The low-recharge seepage rates were observed by seepage meter investigations in this part of the lake (Fig. 11, near Site O2).

Measurements of δ18O in the discharging part of the lake (I2, I3, and offshore piezometers) all indicate inflow. Samples above groundwater values at I3 show that discharge rates at I3 are lower than at I2, which was confirmed by seepage meter data (Fig. 11).

The heterogeneous model shows that a fairly good approximation to direct seepage meter measurement can be made by a large-scale groundwater model, even though seepage meter measurements represent only 0.25 m2 of the lake bed and each model output point represents 2500 m2 of the lake bed. This is especially observed at I2 and I3, where many seepage meter measurements were done. At I3, discharge rates >1 × 10−7 m s−1 were measured at 4 and 6 m of water depth. It is interesting that without these offshore measurements, the simulated discharge from the 50- by 50-m numerical model would have seemed to be overpredicting seepage. Generally, no conclusion can be based on the few offshore measurements made in this study, but they show that if seepage meter rates are to be compared with a regional model, they need to represent the entire lake bed through which the exchanges occur.

Differences in peak discharge and recharge for the homogeneous and heterogeneous models are most probably a result of the reduced interacting area for the heterogeneous model because of the imposed impermeable leakage parameter in the central and western parts of the lake. The higher variance in seepage in the heterogeneous model seems to compare better with the seepage meter measurements.

Ideally the lake bed leakage should be divided into more than the two zones in the heterogeneous model because of the differences in both the sediment composition (gyttja or peat) and the vertical dimensions of the organic layer. The hydraulic conductivity of the organic layer is believed to be below 10−7 m s−1, with strongly decomposed peat values of 10−9 to 10−7 m s−1 (Dahl et al., 2004) and gyttja values in the same range or lower because of its very dense composition. The contrast with the remaining lake bed, with a calibrated value of 6.1 × 10−6 s−1, is so large that changing the organic layer leakage by one order of magnitude in a spatial distributed pattern will not affect the modeled seepage pattern significantly.

Groundwater residence times derived from a groundwater flow model are a result of the simplifications and assumptions associated with the model. One of the simplifications of the Lake Hampen model is that the upper aquifer is only represented by one set of parameters. The hydraulic conductivities for the final model were geometric means representing the whole aquifer. The effect of having homogeneous and steady-state conditions in the upper aquifer is an almost linear relation between depth and simulated mean groundwater residence time. At MLW4, the simulated mean groundwater age below 3.8 m can be described by 1.38 × depth − 1.68 (R2 = 0.997). The expression was calculated from the center depth at Layers 2 to 5 (5.7, 9.5, 13.3, and 19.0 m) and the calculated mean ages from the frequency distribution for the four layers (Fig. 12).

In Fig. 12 it can be seen that CFC ages at the depths of Layers 2 to 5 are fairly close to the simulated ages. The mean simulated ages between 6.6 and 24.9 yr are acceptable compared with the CFC ages. Residence time appears to decrease with depth in Zones 1 and 2, at 2- to 7-m depth. This trend, opposite to what might have been expected, may be real or could be due to the partial degradation of CFC at shallow depths, where dissolved O2 is low (Table 4) and denitrification occurs. For CFC-11, which is most susceptible to degradation (Hinsby et al., 2007), a partial loss of CFC is probably part of the explanation for the high CFC age (Table 4), but for CFC-12, which is a more persistent CFC gas, the observed trend probably reflects reality. The MLW4 well was located on what appeared to be an old lake bottom and observations made during field work indicated lower flow rates in the uppermost layers. Lower hydraulic conductivity in the upper layers may thus be the reason for the observed trend in Zones 1 and 2 (Fig. 12), although mixed flow paths at shallow depths may also play a role.

The reason for the possible low dissolved O2 at shallow levels may be a higher reducing potential compared with deeper levels or it could be due to higher dissolved organic matter in the area where the precipitation originally infiltrated to the aquifer. The reducing potential is no doubt higher in soils of an old lake bottom compared with glacial meltwater sands, which make up the deeper layers, but the extremely high K concentration in the shallow groundwater (Table 4) indicates a large input of organic waste, which is probably responsible for the O2 consumption. The organic waste probably originates from the farm 30 to 60 m upstream and can be traced to a depth of 6 to 7 m. Below Zones 1 and 2, the groundwater chemistry indicates infiltration to the groundwater in agricultural areas, and below Zone 3 (7–11 m) the low NO3 values indicate infiltration from forested areas (Table 4).

Another way of comparing CFC and model ages is to calculate an average CFC age and an average simulated age for the entire observed horizon. If all CFC-12 ages are given equal weights because sampling was done at each meter, the average calculated age of the groundwater reaching the lake (MLW4) is 16 yr. Calculation of an average age for the groundwater model is 13.15 yr when Layer 1 has been given a weight of 0.5 of the other layers (only 50% of the cells in Layer 1 are below the groundwater table and therefore the effective volume is half the size of the other layers).

Matching of measured groundwater and model ages indicate that the model boundary conditions and parameterization simulate the velocity and timing of the aquifer well. This, in consequence, indicates that the discharge of groundwater to the lake is simulated fairly well, which was confirmed by the simulation of seepage rates at the lake bed compared with the seepage meter measurements.

Conclusions and Perspective

Through the combined use of distributed field measurements (hydraulic gradients, seepage meter, GPR, and δ18O), a whole-lake evaluation of recharge and discharge seepage patterns was performed. The study demonstrates that if a reasonable evaluation of the lake–groundwater systems is to be made, whole-lake seepage patterns have to be investigated at different scales. At the near-shore scale, measurements of gradients, and seepage meter rates can be used for detailed studies of local seepage distributions, as shown by Kishel and Gerla (2002), but without results from the offshore parts of the lake, the local-scale results are, at best, difficult to use for managing purposes where the evaluation of lake water or mass budgets are needed. Comparing seepage meter measurements with rates derived from a catchment-scale numerical model is possible if measurements cover the whole range between the lakeshore and the offshore, less permeable areas of the lake. At the intermediate to deeper parts of the lake, investigation of the lake bed by GPR proved to be very effective for mapping low interacting areas at the lake bed (the lake was covered by 32 GPR lines in a day). Knowledge of the lake bed sediments was necessary for conceptualization and modeling of the distributed discharge and recharge patterns at the lake. The addition of the geologic lake bed complexity to the catchment-scale model improved the model performance. The homogeneous model without the lake bed complexity showed that when the geometry of the groundwater-aquifer system is well understood, in this study by a MEP survey, reasonably correct seepage patterns can be simulated by the model.

The recharge zone at Lake Hampen was found by a combination of offshore lake bed mapping, gradient measurements, and identification of the mixing of lake and groundwater by δ18O values observed in deep profiles downgradient of the lake. The recharge pattern seems to be controlled by complex geologic development, with deposition of lake and wetland sediments, in the time since its origin. The study illustrates the importance of knowing the geologic history of a lake to interpret the distribution of seepage and groundwater flow beneath and near the lake. At Lake Hampen, this is exemplified by the important knowledge of the decomposed peat layer formed during a longer period of lower lake water stage, which controls the present recharge pattern.

Finally, the study shows that groundwater ages obtained from CFC studies can be used as data for validation of a groundwater model simulating the interaction of groundwater and surface water. Verification of the modeled groundwater age also validates the modeled groundwater flow rates and thereby the seepage rate from the aquifer to the lake. If the groundwater flow model accurately simulates the age of the groundwater, the seepage rates controlling the water and mass balances of the lake are more likely to be estimated correctly by the model.

We would like to thank Sachin Karan for assistance during the study. The study was supported by the Center for Lake Restoration (CLEAR, Denmark), a Villum Kann Rasmussen Center for Excellence, and the Danish Geocenter (Lake Response).