Estimation of lacustrine groundwater discharge using heat as a tracer and vertical hydraulic gradients – a comparison

Lacustrine groundwater discharge (LGD) can play a major role in water and nutrient balances of lakes. Unfortunately, studies often neglect this input path due to methodological difficulties in the determination. In a previous study we described a method which allows the estimation of LGD and groundwater recharge using hydraulic head data and groundwater net balances based on meteorological data. The aim of this study is to compare these results with discharge rates estimated by inverse modelling of heat transport using temperature profiles measured in lake bed sediments. We were able to show a correlation between the fluxes obtained with the different methods, although the time scales of the methods differ substantially. As a consequence, we conclude that the use of hydraulic head data and meteorologically-based groundwater net balances to estimate LGD is limited to time scales similar to the calibration period.


INTRODUCTION
Lakes are an important part of terrestrial ecosystems and provide essential ecosystem services.Hence, lakes should be in a good ecological state, which is directly related to a proper ecosystem functioning (Baron et al. 2002).Therefore, eutrophication is a major concern in many lakes.One challenge in terms of adequate management is the identification and quantification of the nutrient input paths.Studies focusing on nutrient balances of lakes often neglect input via groundwater, which can hardly be quantified, although some studies showed that lacustrine groundwater discharge (LGD) can play a significant role in the nutrient balance of lakes (e.g.Hayashi andRosenberry 2002, Nakayama andWatanabe 2008).Several authors have spent some effort to determine these fluxes.Sacks et al. (1992) used hydraulic head data to divide the shoreline of the studied lake into areas of exfiltration and infiltration, and used this for numerical modelling.Ommen et al. (2012) used a combination of hydraulic gradients between piezometers and a lake, in combination with seepage meters (Rosenberry and LaBaugh 2008), to assess the nutrient budget of the lake.The use of natural tracers is widely spread and has been successfully applied in many studies.Chemical gradients of natural solutes such as chloride can also be used to estimate fluxes (Hurwitz et al. 2000), as well as temperature profiles in the sediment (Schmidt et al. 2006, Meinikmann et al. 2013).
The method chosen depends on the goal of a study.For studies concerning nutrient balances it is essential to get information about the volume of water entering a lake via LGD and the concentrations of nutrients in this water, as well as the volume of groundwater recharge and nutrient concentrations of the lake water.In a previous study we estimated volumes of LGD and infiltration of lake water into the aquifer for Small Lake Gollinsee and Lake Schulzensee (Rudnick et al. 2014).These calculations were based on a groundwater net balance and vertical hydraulic gradients (VHGs) which were derived from spatial interpolation of piezometer data.The aim of the present study is to compare results from the previous study with discharge rates obtained from inverse modelling of temperature profiles in the lake bed sediments using an analytical solution of the heat transport equation (Bredehoeft and Papaopulos 1965, Schmidt et al. 2006, 2007).

Study site
The lakes which are under investigation in this study are located 50 and 80 km north of Berlin, Germany.Small Lake Gollinsee (N 53.029°, E 13.588°) is located within a former peatland which extends in a north-south direction between ascending terrain in the east and west.The area of the lake is 3.3 ha, the maximum depth is 2.9 m, and the mean depth is 1.8 m.The lake has no surficial inflow or outflow and is a groundwater-fed lake.The lake bed consists of fine organic sediments.The underlying sandy aquifer has a thickness of more than 50 m, according to hydrogeological maps.Lake Schulzensee (N53.247°,E13.275°) is located within a depression in a forested area.The area of the lake is 3.9 ha, the maximum depth is 4.1 m, and the mean depth is 2 m.Like Small Lake Gollinsee, it is a groundwater-fed lake.The lake bed consists of material very similar to that of Small Lake Gollinsee and the underlying aquifer is sandy with a thickness of more than 50 m according to hydrogeological maps.

Groundwater net balance
The groundwater net balance for the studied period from 9 July 2010 to 30 May 2011 (Small Lake Gollinsee) and 12 July 2010 to 30 May 2011 (Lake Schulzensee), was calculated following equation (1): where is the net groundwater input (mm) with positive and negative values indicating net groundwater inflow or outflow, respectively, ΣPrec * is the sum of precipitation (mm), ΣEpot * the sum of potential evaporation (mm), and ∆Lake * the storage change of the lake (mm).Potential Evaporation was assumed to be actual evaporation and was calculated after Allen et al. (1998) for shallow water.The meteorological data needed for this were recorded with automatic weather stations installed on site.Since the stations were not able to detect snowfall for the period 1 December 2010 to 15 March 2011, precipitation data from a station at 20 to 40 km distance were used.In this period the lakes were covered by ice, hence we assumed evaporation to be zero.Changes of lake levels were measured every 4 weeks.

Vertical hydraulic gradients
Based on the assumption that the hydraulic conductivity of the lake bed sediments in comparison to the aquifer is much lower, we calculated VHGs.Monthly records of hydraulic heads of piezometers around the lakes were averaged in the studied time period and spatially interpolated using inverse distance weighting.The spatial distribution of differences between mean hydraulic heads in the aquifer and mean lake level was calculated afterwards.VHGs were derived dividing these differences by the thickness of the sediment layer at the lake bottom (Fig. 1).
Fig. 1 Conceptual drawing of the estimation of spatially distributed VHGs.Hydraulic heads are spatially interpolated between piezometers (h(Aquifer)) and then subtracted by the lake level (h(Lake)).
VHGs are estimated by dividing differences in hydraulic heads (∆h) by the thickness of the lake bed sediments (∆x).

Discharge and recharge rates based on vertical hydraulic gradients
Since steep VHGs can lead to clogging (Rosenberry and Pitlick 2009), we took this into account.
As we could not find any quantitative description of the degree of dependency of clogging on the gradients we introduced a clogging factor as a logistic growth function of the VHG.
Using the following equation ( 2) we derived a spatial distribution of governing hydraulic conductivities: where Kf is the matrix of hydraulic conductivity (m s -1 ), ∆t is the length of the investigated time period (s), Grad is the matrix of VHGs (m m -1 ), Cf is the matrix of clogging factors (-), ΣPrec is the sum of precipitation (m 3 ), ΣEpot is the sum of potential evaporation (m 3 ), and ∆Lake is the storage change of the lake (m 3 ).We used meteorological and averaged VHGs from 9 July 2010 to 30 May 2011 (Small Lake Gollinsee) and 12 July 2010 to 30 May 2011 (Lake Schulzensee) to estimate the spatial distribution of Kf.
To calculate the spatial distribution of LGD and groundwater recharge, we multiplied VHGs by hydraulic conductivities.We did this: (a) with averaged VHGs within the week of temperature profile measurements, and (b) with averaged VHGs over the whole period which was used for finding the Kf distribution.

Discharge and recharge based on temperature depth profiles in the sediment
In July 2010 temperature profiles in the sediment were measured at the lakes.The measurement pattern accommodated the suspected groundwater flow direction.At locations with suspected lacustrine groundwater discharge, five to seven profiles were measured in line, orthogonal to the shore and up to 10 m distance to the shore.An analytical solution of the heat transport equation as shown in equation (3) was used to describe temperature profiles in the sediment for certain discharge rates (Schmidt et al. 2006(Schmidt et al. , 2007)).
where z is the depth of the temperature measurement (m), T(z) is the temperature at depth z (°C) in the sediment, qz is the vertical Darcy-velocity (m s -1 ), pfcf is the volumetric heat capacity of the fluid (J m -3 K -1 ), Kfs is the thermal conductivity of the sediment-water-mixture (J s -1 m -1 K -1 ), L is the depth of the lower boundary condition T(z) = TL (m), TL is the groundwater temperature (°C), and T0 is the temperature of the surface water (°C).We used a parameter set of Kfs = 1.63 J s -1 m -1 K -1 , pfcf = 4190000 J m -3 K -1 , L = 2 m, and TL = 9.7 °C.T0 was chosen individually for each profile.The best fit between measured temperatures in the sediment and calculated temperatures was derived finding a qz,k which minimizes Errork in equation ( 4).
where Errork is the sum of the squared differences between measured temperature Tj in depth zj and calculated temperatures with a Darcy-velocity qz,k.This method is limited to groundwater discharge since it is not possible to estimate a meaningful lower boundary condition for flow of surface water towards the aquifer.Since the locations of the temperature depth profiles were estimated using an ordinary GPS device, the accuracy was low.To compare the results obtained by the use of temperature profiles with the discharge rates estimated based on VHGs, we averaged the results and coordinates of the measurements which were done orthogonal to the shore and adjusted some of the other measurement points by hand.

Vertical hydraulic gradients
Figure 2 shows the spatial distribution of groundwater discharge and recharge in the lakes based on the VHGs within the week of the temperature profile measurements.In the case of Small Lake Gollinsee, groundwater discharges and water infiltration occurred in 53% and 47% of the lake area, respectively.The discharge volume of groundwater is 18 m 3 d -1 and groundwater recharge is -4 m 3 d -1 , which adds to a net groundwater inflow of 14 m 3 d -1 .The net groundwater balance based on meteorological data for the period from 6 to 9 July 2010 yields a net inflow of groundwater of 49 m 3 d -1 .We take the meteorological-based groundwater net balance as a reference since we used this to find the optimum values for the hydraulic conductivities of the lake bed sediments.Hence, the fluxes, which we obtained for the week.underestimate the inflow of groundwater by approx.72%.
In the case of Lake Schulzensee some piezometers at the eastern shore of the lake showed high hydraulic heads which lead to strong gradients and therefore high fluxes.Of the lake area, 91% are considered to be zones of groundwater discharge and 9% to be recharge areas.The groundwater discharge volume of 247 m 3 d -1 and only -6 m 3 d -1 surface water infiltration leads to a net groundwater balance of 241 m 3 d -1 .This result is contradictory to the meteorological-based groundwater net balance, which suggests a net loss of water to the aquifer of 50 m 3 d -1 .This is mainly a time-scale issue, since we calibrated the hydraulic conductivity of the lake bed sediments to data which was averaged over a period of roughly one year.

Temperature depth profiles
The adjusted locations of measured temperature profiles in the lakebed sediments are shown in Fig. 2. The groundwater discharge rates at these locations range from 9 to 27 L m -2 d -1 with a mean of 20 L m -2 d -1 (Small Lake Gollinsee) and 11 to 42 L m -2 d -1 with a mean of 23 L m -2 d -1 (Lake Schulzensee).In the case of Lake Schulzensee three locations which suggest surface water infiltration could be identified near the southwestern shore.Since quantification of these fluxes is not possible we neglected them.We could not find temperature profiles which indicate surface water infiltration at Small Lake Gollinsee, although differences between hydraulic heads of the piezometer and lake level clearly suggested influent conditions in the northwestern part of the lake.This could be caused by strong non-stationary conditions that affect the method we used.

Comparison of results from the two methods
A comparison of the results obtained by the two different methods described above is presented in Fig. 3.For each temperature profile location the corresponding flux estimated using VHGs was determined.We used two different flux calculations: (a) mean VHGs of the week when temperature profiles were measured (indicated with circles in Fig. 3), and (b) mean VHGs over the whole studied period of roughly one year (indicated with triangles in Fig. 3).In general, the fluxes which were estimated using temperature profiles are higher than the values based on VHGs.Nevertheless, a fair correlation between the fluxes determined by the two methods can be found with r = 0.54 (week) and r = 0.67 (year) for Small Lake Gollinsee and r = 0.41 (week) and r = 0.67 (year) for Lake Schulzensee.The higher correlation of temperature profile based values with the mean annual fluxes could be led back to the fact that calibration was done on an annual basis.In the case of Small Lake Gollinsee the year-based fluxes are higher than the week-based; in the case of Lake Schulzensee it is the opposite.Some of the variation of the temperature profile-based fluxes can be explained by the high spatial heterogeneity which is often observed when dealing with surface water-groundwater interaction.This heterogeneity cannot be reproduced in the VHG approach.Furthermore, the use of the analytical solution of the heat transport equation can lead to a significant overestimation of groundwater discharge when the assumption of stationarity is violated (Anibas et al. 2009, Schornberg et al. 2010).

CONCLUSION
Although the results obtained by these two methods differ, the patterns concur, which is shown by statistical correlation.Limitations of the approach to estimate groundwater discharge and recharge using interpolated VHGs are the difficulty to reproduce spatial heterogeneity and the dependency on the temporal scale of calibration procedures.This means that the use of hydraulic heads for the estimation of groundwater discharge and recharge is not appropriate to assess short-term balances, but gives reasonable results on the time scale of calibration periods.

Fig. 2
Fig. 2 Fluxes calculated using mean VHGs in the week of the temperature measurements.Contour lines indicate levels of fluxes.Squares and asterisks indicate locations of the piezometers and sediment temperature profiles respectively.

Fig. 3
Fig.3Fluxes based on VHGs plotted versus fluxes obtained by using heat as a tracer.Plotted with circles: VHGs averaged during the week where temperature profiles were measured, triangles indicate fluxes from hydraulic heads averaged over the study period.The dashed line indicates a perfect agreement.The correlation coefficient is for Small Lake Gollinsee r = 0.54 (week) and r = 0.67 (year); for Lake Schulzensee r = 0.41 (week) and r = 0.67 (year).