Quantifying Groundwater Dependence of a Sub-polar Lake Cluster in Finland Using an Isotope Mass Balance Approach

A stable isotope study of 67 kettle lakes and ponds situated on an esker aquifer (90 km 2) in northern Finland was carried out to determine the role and extent of groundwater inflow in groundwater-dependent lakes. Distinct seasonal fluctuations in the δ 18 O and δ 2 H values of lakes are the result of seasonal ice cover prohibiting evaporation during the winter. An iterative isotope mass balance approach was used to calculate the inflow-to-evaporation ratios (I TOT /E) of all 67 lakes during the summer of 2013 when the isotopic compositions of the lakes were approaching a steady-state. The balance calculations were carried out independently for 2 H and 18 O data. Since evaporation rates were derived independently of any mass balance considerations, it was possible to determine the total inflow (I TOT) and mean turnover time (MTT) of the lakes. Furthermore, the groundwater seepage rates to all studied lakes were calculated. A quantitative measure was introduced for the dependence of a lake on groundwater (G index) that is defined as the percentage contribution of groundwater inflow to the total inflow of water to the given lake. The G index values of the lakes studied ranged from ca. 39 to 98 %, revealing generally large groundwater dependency among the studied lakes. This study shows the effectiveness of applying an isotope mass balance approach to quantify the groundwater reliance of lakes situated in a relatively small area with similar climatic conditions.


Introduction
The characterization of groundwater dependent ecosystems (GDEs) is a requirement of the Groundwater Directive (EC, 2006).These systems are often complex and their hydrology and contact with aquifers are not well established.Lakes can be dependent on groundwater directly or indirectly, and this dependence can vary over time (Kløve et al., 2011).Understanding groundwater and lake water interaction is important not only for water resource management (Showstack, 2004), but also for understanding the ecology and eutrophication of lakes, since groundwater may be a key element in the lake nutrient balance (Ala-aho et al., 2013;Belanger et al., 1985;Brock et al., 1982;Kidmose et al., 2013).Furthermore, the vulnerability of lakes to pollution can be controlled by their dependency on groundwater (Kløve et al., 2011).Methods such as seepage meters (Ala-aho et al., 2013;Rosenberry et al., 2008Rosenberry et al., , 2013)), environmental tracers (e.g.Dinçer, 1968;Shaw et al., 2013;Stets et al., 2010;Yehdegho et al., 1997;Zuber, 1983) and numerical modelling (e.g.Krabbenhoft et al., 1990;Stichler et al., 2008;Winter and Carr, 1980) can be used to determine the groundwater reliance of lakes.
Heavy stable isotopes of water ( 18 O, 2 H) can be considered as ideal tracers for studying the hydrological cycle (e.g.Clark and Fritz, 1997).Fractionation of isotopes of water is the very factor enabling their use in hydrological studies, as it governs the changes in isotopic abundances within the water cycle (Gat, 2010).At a global scale, the 2 H and 18 O isotope composition of meteoric waters cluster along the line called the global meteoric water line (GMWL), first determined by Published by Copernicus Publications on behalf of the European Geosciences Union.E. Isokangas et al.: Quantifying groundwater dependence of a sub-polar lake cluster Craig (1961): δ 2 H = 8 • δ 18 O + 10.Locally, this linear relationship may have a slightly different form (local meteoric water line -LMWL).Evaporation from an open water body fractionates isotopes so that the remaining liquid phase is enriched in both 2 H and 18 O in proportion with their effective fractionation factors accompanying this process.Consequently, the isotopic composition of the evaporating water body evolves in the δ-space along the line known as the local evaporation line (LEL), whose slope is significantly smaller than that characterizing the local or global meteoric water lines.The position of the isotopic composition of lake water along this line is strongly related to the water balance of the lake (e.g.Gat, 1996;Gibson and Edwards, 2002;Rozanski et al., 2001).
The methodology of isotope-aided studies of the water balance of lakes has been thoroughly discussed in a number of review papers and textbooks (e.g.Darling et al., 2005;Froehlich et al., 2005;Gat and Bowser, 1991;Gat, 1995;Gonfiantini, 1986;Rozanski et al., 2001).Although several authors have applied isotope techniques in studying lakes in cold climates (Gibson and Edwards, 2002;Gibson, 2002;Gibson et al., 1993;Jonsson et al., 2009;Turner et al., 2010;Yi et al., 2008), mostly in Canada and northern Sweden, these studies were generally focused on lakes spread over large areas.
The central aim of this study was to quantify the groundwater dependence of 67 kettle lakes and ponds situated across a relatively small area (90 km 2 ) of the Rokua esker aquifer occupying a large glaciofluvial deposit in northern Finland.To quantify the extent of the interaction between the aquifer and the lakes, a dedicated isotope study was launched in 2013.This was part of comprehensive investigations (2010)(2011)(2012) aimed at understanding the hydrology of an esker aquifer area where some of the kettle lakes and ponds have suffered from water level decline or eutrophication.Since the seasonal isotopic behaviour of the selected lakes in the study area was already fairly well understood based on the data collected from 2010 to 2012, it was decided to conduct a large-scale one-time survey of the isotopic composition of all 67 lakes on the Rokua esker in order to quantify their dependence on groundwater.Ala-aho et al. (2013), who studied 11 lakes on the esker, showed that the water levels in closedbasin seepage lakes have more fluctuations than the drainage lakes, which have more stable water levels.On the other hand, the drainage lakes are more productive.Their study also showed that subsurface flow can transport phosphate to lakes.Therefore, it was important to quantify the groundwater dependence of all lakes on the esker and propose an appropriate index reflecting this dependence.The large-scale field campaign conducted in July and August 2013 comprised the sampling of water in all 67 lakes for isotope analyses, combined with continuous temperature measurements and aerial thermal imaging of the lakes.

The study area
The Rokua esker aquifer area, situated in northern Finland, was formed during the transition period from Late Glacial to Holocene, between approximately 12 000 and 9000 years ago (Tikkanen, 2002).As ice retreated, a long ridge formation, consisting mainly of fine and medium sand (Pajunen, 1995), was shaped.Ancient sea banks surrounding the esker show that the esker was originally an island that gradually rose from the sea (Aartolahti, 1973).Today, the highest elevation of the esker is 100 m above the surrounding low-lying peatlands and the layer thickness of sand ranges from 30 m to more than 100 m above the bedrock.Sea banks, dunes and kettle holes form a rolling and geologically unique terrain (Aartolahti, 1973).Kettle holes were formed when ice blocks were buried in the ground and, as they melted, left depressions in the landscape.The ground surface of the esker is mainly lichen-covered pine forests.Hydrologically, the Rokua esker is an unconfined aquifer, one of the largest in Finland, and it has two regional groundwater mounds (Rossi et al., 2014).The recharge area of the aquifer is 90 km 2 and the discharge zones are situated in the surrounding peatlands, which partially confine the aquifer (Rossi et al., 2012).
Kettle holes -long and narrow depressions -give Rokua esker its distinct character.The sizes of these kettle holes vary.They can be 1-80 m deep, between 10 m and 1.5 km long, and 0.4 km wide (Aartolahti, 1973).Most of the kettle holes are now dry, but due to the influence of groundwater in the past, peat has accumulated at the bottom of them, creating kettle hole mires (Pajunen, 1995).However, the alternating topography of the area is reflected in the existence of approximately 90 lakes or ponds, referred to as kettle lakes or ponds.Peat started to accumulate in the border regions of the lakes more than 8000 years ago, so most of the kettle lakes and ponds are partly paludified (Pajunen, 1995).Nevertheless, the majority of the lakes and ponds are characterized by their crystal clear water, which attract people; a number of holiday homes and hotels are located on the lakeshores.The lakes are nowadays widely used for different recreational activities, such as swimming, fishing and scuba diving (Anttila and Heikkinen, 2007).The uniqueness of the glaciofluvial formation of Rokua, in which the actions of ice, water and wind can be seen, has been recognized in many ways.Some of the Rokua esker is protected by Natura 2000 and by the Finnish nature reserve network.Rokua was recently chosen to be part of the UNESCO GeoPark Network and is currently the northernmost region in this network.yse the stable isotopic composition of water, nutrients, water quality parameters (T , pH, EC, O 2 ) and geochemical parameters (silica, major cations and anions) (Fig. 1).During the field campaign conducted in July and August 2013, a total of 67 lakes and ponds were surveyed for the same parameters and thermal images of lakes taken from the air in a helicopter using a FLIR thermal camera.In addition, composite (monthly) precipitation samples were collected during the open water season at the station on the esker during 2010-2013.Precipitation samples for winter were collected once a year before snowmelt by taking a uniform sample of the whole snowpack depth.Water quality parameters were analysed in the field using WTW Multi 3430 or Multi 350i meters for temperature, oxygen, EC and pH.Samples of lake water were collected with a Limnos sampler, approximately 1 m below the water surface and 1 m above the bottom of the lake.If the depth of the lake was less than 2 m, only one sample from the depth of 1 m was taken and if it was more than 20 m, samples were taken from the middle of the water profile as well.Depending on their shape and size, the lakes had between 1 and 4 sampling locations.Stream samples were collected by submerging a bottle in water, facing upstream.Piezometers were pumped for at least 10 min prior to taking groundwater samples or until the colour of the water was clear.The samples were collected one metre below the water table.All sampling bottles (HDPE) were rinsed with the sampled water prior to filling.

Hydrological measurements and thermal imaging
Samples for isotope analyses were stored in the dark at a reduced temperature (4 The isotopic composition of water samples was analysed using CRDS technology with a Picarro L2120-i analyser.Samples with visible colour or suspended matter were filtered (pore size 25 µm) prior to analysis.The measured 2 H / 1 H and 18 O / 16 O isotope ratios are reported as relative deviations from the VSMOW standard.Typical uncertainty of the reported δ 18 O and δ 2 H values are ±0.1 and ±1.0 ‰ respectively.
Lake water temperature was measured continuously during the ice-free period in 2013 for two lakes on the Rokua esker: Ahveroinen (3.3 ha) and Saarinen (15.3 ha).Hobo loggers (pendant temperature data logger UA-001-08 and conductivity logger U24-001, accuracy 0.1 • C) were installed 50 cm below the lake surface.In addition, the surface water temperature for lake Oulujärvi (92 800 ha) was obtained from the database of the Finnish Environmental Institute (2013).Lake Oulujärvi is located east, next to the study site, 1 km from the easternmost lake studied.Thermal imaging of the lakes was conducted on 5 August 2013 using a Flir Thermacam P-60 thermal camera.This camera had 320 × 240 pixel sensor resolution and an opening of 24 • .It covered the electromagnetic spectrum from 7.5 to 13 µm.The imaging was taken by helicopter 150 m above the lakes.The image data were correlated to the predominant weather conditions (temperature and relative humidity) with data from the FMI Pelso weather station measured every 10 min.

E. Isokangas et al.: Quantifying groundwater dependence of a sub-polar lake cluster
Depth profiling was undertaken in the lakes for which no depth contour lines are available (National Land Survey of Finland, Basic map, https://tiedostopalvelu. maanmittauslaitos.fi/tp/kartta).It was carried out with a portable depth-sounding radar (resolution 0.1 m) or with a measuring cable and GPS system.Typically, two profiles from the shore to the deepest point were defined.The number of measurement points differed between the lakes depending on their size.In total, 52 lakes were surveyed.

Lake volumes
The volumes of the lakes were determined in an ArcGIS environment using depth profiling measurements, contour lines and border lines.The water surface levels of the lakes were estimated using elevation levels presented in the basic map by the National Land Survey of Finland.Lake morphology was mostly interpolated using a spline that results in a smooth surface passing by all the input points (ESRI, 2014).The tension method with 0.1 weight and 3 input points was used to calculate the values for the interpolated cells.Interpolation rasters were extracted by surface water areas.The mean depths of these new rasters were multiplied by the water surface areas in order to calculate the volumes of the lakes.

Evaporation from lakes
Evaporation (E) was calculated individually for all the lakes surveyed using a mass transfer approach (Rosenberry et al., 2007;Dingman, 2008).This method was chosen because it yields instantaneous rates of evaporation and takes into account lake sizes (Harbeck, 1962).Our data enabled the calculations of daily mean values of evaporation flux for all the lakes studied.Averaging the parameters for periods of time longer than one day can lead to significant biases in the calculated E values since on these time scales the vapour pressure differences and wind speed may correlate (Jobson, 1972).The following expression was used to calculate the evaporation flux (Dingman, 2008): where -E is evaporation rate (mm s −1 ); -K E is the mass-transfer coefficient in m km −1 kPa −1 describing the impact of turbulent eddies of the wind on vertical transport of water vapour from the lake with area A L (km 2 ), K E = 1.69 × 10 −5 • A −0,05 L (Harbeck, 1962); v a is wind speed (m s −1 ) at 2 m height; e s is the saturation vapour pressure in kPa at surface water temperature T s ( • C), e s = 0.611 • exp 17.3•T s T s +237.3 ; e a is the vapour pressure in the air in kPa, e a = h•0.611•exp 17.3•T a T a +237.3 , where h is relative humidity and T a is air temperature ( • C).
Wind speed measured at 10 m height was adjusted to the corresponding speed at 2 m height using the power law profile (Justus and Mikhail, 1976): where v r is the measured wind speed at the reference height z r (10 m), z is the height for which speed is adjusted (2 m) and β is the friction coefficient.The value of 0.15 for β, characteristic for grassland, was employed in our study (0.1 characterizes oceans and lakes) (Bañuelos-Ruedas et al., 2010) since the lakes are surrounded by forests lowering the wind speed.
The meteorological parameters necessary for the calculations (relative humidity, wind speed and air temperature) were obtained from the meteorological station 5502 (Vaala-Pelso) of the Finnish Meteorological Institute (2014), located approximately 10 km from the site.A probable range for the lakes' surface water temperature was evaluated using continuous temperature measurements at 50 cm depth from one of the studied lakes, Ahveroinen (see Sect. 4.1), and a standard deviation of this temperature determined from thermal images.Using the derived temperature range, a probable range for evaporation rates from all the lakes was calculated.The adopted method relies only on temperature difference measured on one day, but since all the lakes are in a relatively small area with almost identical weather conditions, it is highly probable that the seasonal behaviour of the surface water temperature of the studied lakes is similar.

Isotope mass balance
Instantaneous water and isotope balances for an evaporating surface water body can be formulated as follows: where V stands for the volume of the surface water body, δ L signifies its isotopic composition and I TOT , E and O TOT represent the total inflow, evaporation and total outflow of water from the system, respectively, whereas δ IT , δ E and δ OT stand for their respective isotopic compositions, expressed in ‰.
As the total inflow may consist of several components (precipitation, underground and surface inflows), each with its specific isotopic composition, δ IT should be calculated as a flux-weighted mean of the respective isotopic compositions of individual components.The total outflow may also consist of surface and underground components.For well-mixed systems it is typically assumed that δ OT = δ L .The isotopic composition of the evaporation flux, δ E , cannot be measured directly.However, it can be calculated using the expression derived from the linear resistance model describing isotope effects accompanying the evaporation process (Craig and Gordon, 1965;Horita et al., 2008): where ε is the total effective isotope fractionation, ε = ε * + ε, where ε * stands for equilibrium isotope enrichment (ε * = (1 − 1/α LV ) × 10 3 ) expressed in ‰ and α LV represents the equilibrium isotope fractionation factor between liquid and gaseous phases.The kinetic isotope enrichment ε, is defined as where C k stands for the kinetic enrichment parameter.δ A represents the isotopic composition of atmospheric water vapour over the evaporating water body (‰) and h N is the relative humidity of the local atmosphere, normalized to the temperature of evaporating water.
When the evaporating water body is in hydrologic and isotopic steady-state (dV / dt = 0 and dδ L / dt = 0, respectively), the following approximate expression describing the isotope enrichment of the evaporating water body can be derived from Eqs. (3-5) (see e.g.Gat and Bowser, 1991): where δ stands for the evaporative enrichment and δ LS is the steady-state isotopic composition of the studied system.
The following expression describing the ratio of the total inflow to the evaporation rate can be derived from Eq. ( 6): If the evaporation flux E can be assessed independently, as is the case of this study, the total inflow of water to the given lake can be calculated from Eq. ( 7).From this, the groundwater component of this inflow can be further inferred, provided that the other components of this inflow are known or can be independently assessed.
To quantify the water balance of a lake with the aid of Eq. ( 7) the knowledge of the isotopic composition of atmospheric moisture interacting with the lake, δ A , is required.Field measurements of this parameter have become feasible only recently thanks to advancements in CRDS technology.However, the combined studies of the isotopic composition of atmospheric water vapour and precipitation performed in moderate climates (Jacob and Sonntag, 1991;Schoch-Fischer et al., 1984) have shown that, on a monthly basis, the isotopic composition of precipitation is generally in isotopic equilibrium with the local atmospheric moisture at groundlevel temperature.This is particularly true for the summer season.Therefore, the value of δ A can be derived from the isotopic composition of local precipitation, δ P , which is also required to quantify the isotopic composition of the total inflow to the studied lake, δ IT , appearing in Eq. ( 7).The following relation can be used to calculate δ A (in ‰): Equation ( 7) is valid under two basic assumptions: (i) the evaporating water body is in hydrologic and isotopic steadystate, and (ii) the water body is isotopically homogeneous.Natural surface water systems, such as lakes, typically operate close to their hydrologic and isotopic steady-states attained in the course of their long history.Their steady-state characteristics are defined by local climate, morphological setting and prevailing hydrological regime.Such systems usually exhibit seasonal fluctuations of varying amplitude, caused by seasonal fluctuations of local climate (surface air temperature, relative humidity, precipitation amount), superimposed on long-term trends.Seasonal ice-cover, typical for mid and high latitudes, may also contribute to the seasonal fluctuations of the steady-state characteristics of such systems.The gradual attainment of the isotopic steady-state, which is characterized by an exponential function describing the temporal evolution of δ L , can be observed only for surface water systems which are artificially created, such as dredging lakes resulting from exploitation of gravel deposits (Zimmerman, 1979).The time constant characterizing the dynamics of this process is mainly controlled by the mean turnover time of water in this system, defined as the ratio of its volume to the total inflow, its hydrological balance (I TOT /E ratio) and the normalized relative humidity (e.g.Gonfiantini, 1986;Zimmerman, 1979).Lake Oulujärvi Lake Ahveroinen 1 Lake Saarinen Surface air temperature Figure 3. Daily mean surface water temperatures of lakes Oulujärvi (92 800 ha) (Finnish Environmental Institute, 2013), Ahveroinen 1 (3.3 ha) and Saarinen (15.3 ha) during the summer of 2013, compared with the surface air temperature data for the same period.Lake Oulujärvi is located in the east, next to the study site, 1 km from the easternmost lake studied.Vertical lines mark the period used in the calculations of evaporation and isotope mass balance.

Climate and lake water temperature data
The climate is strongly seasonal in the Rokua esker study area.The long-term monthly mean values of surface air temperature vary from −10.9 • C (January) to 13.2 • C (June) for the period between 1959 and 2013 (Fig. 2).The amount of monthly precipitation varies from 29 mm (February and April) to 79 mm (August) for the same period.The warmest months of the year are June, July and August.The long-term  monthly mean relative humidity of air varies from 61 % (May) to 91 % (November).
The seasonal temperature patterns of the monitored lakes were very similar, despite significant differences in lake size (Fig. 3).Thus the surface water temperature of lake Ahveroinen 1 (mean value of 19.1 • C) during the period from 1 June 2013 to 31 August 2013 was used as a basis for estimating the water temperature of other lakes.Thermal images collected on 5 August 2013 yielded comparable surface water temperatures that ranged from 19.5 to 24.6 • C, with a mean of 21.3 • C and a standard deviation of 0.87 • C. Combining the results of continuous temperature measurements and thermal images, an estimate of the mean surface water temperature of all lakes for the period from 1 June to 31 August 2013 was derived to be 19.10 • C ±0.87 • C.This temperature was used in the isotope mass balance calculations.

Local isotopic compositions
An overview of the isotopic composition of different types of water in the study area is presented on the δ 2 H-δ 18 O space in tion located on the esker, the mean isotopic compositions of the selected lakes, streams and groundwater monitored during the period 2010-2013 ( Although majority of groundwater samples cluster near the LMWL-LEL intersect, there are some data points lying along the LEL in the δ 18 O-δ 2 H plot indicating the contribution of (evaporated) lake water to groundwater.However, the interconnection between the lakes via groundwater a calculated for the period 1 June-31 August 2013; b lakes with identified surface inflow from an upstream lake; c since the upstream lake for lake Nimisjärvi was not sampled, the mean isotopic composition of total inflows to lakes with identified surface water inflows was used for isotope mass balance calculations of this particular lake; d as the isotopic composition of a lake becomes comparable with the isotopic composition of the total inflow, the isotope mass balance calculations become very uncertain.Therefore, for lake Siirasjärvi 2 ( δ 18 O = 0.45 ‰ and δ 2 H = 1.6‰) the values of I TOT /E, MTT and G index reported in the table have only an indicative character.is likely to be minor since lakes probably have groundwater table maxima between them as they are situated in deep holes in the landscape.Nevertheless, the existence of deeper flow paths from the upper elevation lakes to the lower ones cannot be excluded.Based on the estimation of elevation differences and distances between the lakes, we identified the lakes which do not have surface inflows but which might receive some groundwater input originating from upper elevation lakes.These are lakes No. 6,14,20,34,38,41,47,50 and 59.The sensitivity study presented in Sect.4.5 considers probable changes of the isotopic composition of the total inflow to each lake caused by the presence of evaporated lake water component.The influence of evaporated lake water seeping to groundwater is illustrated in Fig. 5, showing the isotopic composition of lake Ahveroinen 1 and adjacent groundwater.The mean δ 18 O and δ 2 H values of piezometers MEA 2010 and MEA 1907 situated on the south-eastern and northwestern sides of the lake were −13.4 and −97 ‰, and −10.7 and −83 ‰ respectively, clearly indicating a substantial (ca.55 %) contribution of lake water to groundwater at the northwestern side of the lake.A smaller contribution (ca. 10 %) of lake water to groundwater can be seen on the eastern and western sides of the lake where the mean δ 18 O and δ 2 H values of groundwater were −12.9 and −93 ‰ respectively.The main direction of groundwater flow is therefore from southeast to north-west, which coincides with the results from seepage measurements conducted by Ala-aho et al. (2013).The direction of groundwater flow can also be noted from the difference in the mean isotopic composition of the lake water between points 2 and 3 of −8.7 and −73 ‰, and −8.5 and −72 ‰ respectively.The difference in isotopic compositions between points 2 and 3 was greatest during the winter: in March 2011 the difference in δ 18 O and δ 2 H between these points was −1.1 and −4 ‰ respectively.

Temporal variations in the isotopic composition of lake water
The seasonal variability in the isotopic composition of the lakes studied is illustrated in Fig. 6, showing changes of δ 18 O in lake Ahveroinen 1 at two depths (1 and 4 m).δ 18 O of lake Ahveroinen 1 reveals distinct seasonal fluctuations with peak-to-peak amplitude in the order of 1 ‰.This lake does not have any surface inflows or outflows.After the disappearance of ice cover during the spring (April-May), the lake starts to evaporate, which results in its gradual enrichment in heavy isotopes, approaching the steady-state value sometime in September-October.Freezing of the lake in late autumn stops the evaporation flux.Systematic decline of δ 18 O during the ice-cover period seen in Fig. 6 stems from gradual dilution of lake water with groundwater seeping into the lake.
Figure 6 shows that the lake is well mixed throughout the year.
The declining parts of the δ 18 O curve in Fig. 6 can be used to assess the intensity of groundwater inflow during the icecover period, if the volume of the studied lake is known and the isotopic composition of groundwater is constant.The isotope balance of such a lake system (Eq.4) can be then expressed as follows: Since I TOT = I GW , δ IT = δ GW , O TOT = O GW = I GW and δ OT = δ L , Eq. ( 9) becomes: The solution of this differential equation reads as follows: where k = I GW /V and δ Lo is the isotopic composition of the lake at the beginning of the ice-cover period (δ 18 O Lo = −8.5 ‰ for lake Ahveroinen 1).Equation ( 10), with δ L = δ Lo , can be used to calculate the mean flux of groundwater to the lake during the ice-cover period.The observed reduction of δ 18 O of lake Ahveroinen 1 by ca.1.2 ‰ over the six-month period (see Fig. 6, Table 1) results from the continuous inflow of groundwater with specific isotopic signature (δ 18 O GW = −13.4‰).The groundwater seepage rate obtained from Eq. ( 12) is ca.160 m 3 day −1 .Identical value was obtained using corresponding 2 H data.It is worth noting that groundwater inflow to lake Ahveroinen 1, derived for the summer period of 2013 from isotope mass balance calculations (ca.300 m 3 day −1 ) is almost two times higher, which suggests significant seasonal variations of groundwater inflow to Rokua lakes, with high groundwater fluxes during summer and low fluxes during winter.These variations most probably stem from seasonal fluctuations in the water table caused mainly by snow melt induced groundwater recharge in spring.Observations for lake water stage and adjacent piezometer MEA 2010 (see Fig. 5) during two years demonstrated an average 20 % higher hydraulic gradient in the summer than in the winter.This is in agreement with our data showing higher groundwater inflow during summer.The annual variability in groundwater fluxes to lakes in the study area was also shown in a modelling study by Ala-aho et al. (2015).

Quantifying groundwater dependence of the studied lakes
Out of the 67 Rokua lakes sampled during the July-August 2013 field campaign, 50 lakes did not reveal any surface water inflow in the form of a stream or creek.For the remaining 17 lakes, surface inflows were identified.For all but one lake those surface inflows could be linked to specific upstream lakes which were sampled during the July-August 2013 campaign.
The isotopic composition of the lakes sampled covers a wide range of δ values: from −5.6 to −12.7 ‰ and from −57 to −93 ‰ for δ 18 O and δ 2 H respectively.This large variability reflects a wide spectrum of the heavy isotope enrichment of the lakes studied.Since the Rokua lakes are situated in a unique climatic region, their observed isotopic composition is primarily controlled by their water balance, which in turn can be characterized by the total inflow-to-evaporation ratio.Equation ( 7) was used to calculate the total inflow-toevaporation ratios (I TOT /E) for all studied lakes.The isotope mass balance calculations were run for the period of June-August 2013 separately for 18 O and 2 H. Since the evaporation rates from the studied lakes were derived from Eq. ( 1), independently of any mass balance considerations, the total inflow to each lake could also be calculated from the assessed I TOT /E ratios.
The values of the parameters occurring in Eq. ( 7) were derived as follows.
i.The mean isotopic composition of atmospheric moisture (δ A ) was calculated with the Eq. ( 8) using the available isotope data for local precipitation.The following mean δ A values were used in isotope mass balance calculations: δ 18 O A = −20.4‰ and δ 2 H A = −149 ‰.
ii.The mean relative humidity, normalized to lake water temperature (h N ) was calculated using the mean surface air temperature (+14.9 • C), the mean surface water temperature of the lakes (+19.1 • C) and the mean relative air-based humidity calculated on the basis of daily mean values available from meteorological station in the area (79.2 %).The resulting mean h N value was 60.7 %.
iii.The total effective isotope fractionation (ε) was derived as the sum of the equilibrium and kinetic isotope enrichments (ε = ε * + ε).The equilibrium isotope enrichment (ε * ) values for 18 O and 2 H, were calculated for the mean surface water temperature of +19.1 • C using the known temperature dependence of the empirical equilibrium fractionation factors, α LV (Horita and Wesolowski, 1994).The values of kinetic enrichment parameters (C k ) used to calculate ε, 14.2 ‰ for 18 O and 12.5 ‰ for 2 H, were adopted after Gonfiantini (1986).Those values were obtained in wind-tunnel experiments (Vogt, 1976) and are widely used in lake studies.iv.Isotopic homogeneity of the studied lakes was addressed through multiple samplings of large lakes (in both the horizontal and vertical direction -see Sect.3.1).The range of the measured δ 18 O values for the given lake was generally lower than one per mil.It was assumed that the average values calculated on the basis of individual measurements performed in each lake represent the studied systems sufficiently well.
v. The isotopic composition of lake water obtained during the sampling campaign in July and August 2013 for each studied lake was used in the isotope mass balance calculations.As discussed in Sect.3.4 above, the isotopic compositions of the studied lakes fluctuate seasonally reaching their steady-state values toward the end of the ice-free period (September-October).Therefore, the isotopic compositions of lake water samples collected during the late summer (July-August) may still deviate slightly from the respective steady-state values, thus creating uncertainty in the assessed components of the water (see Sect. 4.5).
vi.As the isotopic compositions of the surface and underground components of water inflow to each studied lake were not measured directly, an iterative approach was adopted to calculate δ IT individually for each lake.In the first step it was assumed that δ IT is defined by the intercept of the LMWL and LEL lines (Fig. 4).With these δ IT values the underground component of I TOT for each lake was derived from Eq. ( 7).The second step differed for the lakes without surface inflow and with identified surface inflow from an upstream lake.For the lakes without surface inflow, the δ IT values were calculated individually as flux-weighted averages of the underground inflow obtained in the first/previous step of the procedure and precipitation input to the given lake, each with their respective mean isotopic compositions representing the period of June-August 2013 (see Sect. 4.2 above).For the lakes with identified surface inflow from an upstream lake, in the first instance the total outflow of water from the upstream lake was calculated using the appropriate mass balance equation.Then, it was assumed that 25 % of the total outflow from the upstream lake flows to the downstream lake as surface inflow carrying the isotopic composition characteristic for the upstream lake.The δ IT values were calculated individually for each lake belonging to this group of lakes as flux-weighted averages of three components: precipitation, surface inflow and groundwater inflow.For both lake groups, the calculations were repeated until the change of δ IT in subsequent iteration step was in the order of the analytical uncertainty of isotope measurements (0.1 ‰ for δ 18 O and 1 ‰ for δ 2 H).The calculated inflow-to-evaporation ratios of the studied lakes based on 18 O isotope mass balance are shown in Fig. 7 as a function of isotope enrichment of lake water with respect to the isotopic composition of the total inflow to the given lake.They vary in a wide range, from I TOT /E values between 2 and 3, and large 18 O isotope enrichments between approximately 6.5 and 8.0 ‰, indicating evaporation dominated systems, to typical through-flow lakes characterized by I TOT /E ratios higher than 10 and moderate 18 O isotope enrichments of less than 2 ‰ (Table 2).Knowing the volume and the total inflow, the mean turnover time of water (MTT) in each lake could be quantified as the ratio of lake volume to the total inflow.The calculated MTT values range from approximately one week for the Pasko pond (V = 2 × 10 3 m 3 , mean depth 0.2 m, maximum depth 0.7 m) to approximately five years for lake Saarijärvi 2 (V = 2.47 × 10 6 m 3 , mean depth 11.8 m, maximum depth 26 m), with the mean in the order of ten months (Table 2).Lake Saarijärvi 2 is the deepest of all the lakes surveyed.As expected, the calculated MTT values correlate well with the mean depth of the studied lakes, expressed as the volume-to-surface area ratio.However, the link between MTT and the I TOT /E ratio is much weaker; lakes with higher I TOT /E ratios tend to have shorter mean turnover times.
The dependence of the studied lakes on groundwater can be quantified through an index (G index) defined as the percentage contribution of groundwater inflow to the total inflow of water to the given lake.Groundwater inflow was derived by subtracting the precipitation and surface water inflow (if it exists) from the total inflow.Such an evaluation was undertaken for all the lakes listed in Table 2.Note that for the group of lakes with identified surface inflow from an upstream lake it was assumed arbitrarily that this surface inflow is 25 % of the total outflow from the upstream lake (discharges of surface inflows were not measured).The resulting groundwater seepage rates vary from less than 20 m 3 day −1 for Kissalampi pond to around 14 × 10 3 m 3 day −1 for lake Nimisjärvi, the lake with the largest surface area (167.5 ha) among all the studied lakes.
Figure 8 summarizes the values of G index obtained for the lakes surveyed during the July-August 2013 sampling campaign.The mean values of G index obtained from 2 H and 18 O balance are shown.They vary from ca. 40 to more than 95 %.The lowest value (39.4 %) was obtained for lake Etu-Salminen.The highest G values were derived for lakes Kiiskeroinen (97.1 %) and Levä-Soppinen (97.5 %).Interestingly, these lakes are characterized by a high degree of eutrophication induced by high loads of phosphorus brought to the lakes with groundwater (Ala-aho et al., 2013).Although the G index describes the groundwater dependency of the studied lakes rather unambiguously, lakes with moderate G index values could also suffer if the groundwater table in the esker aquifer were to decline as a result of climate and/or land-use changes.
The isotope mass balance calculations for Rokua lakes were run independently for 18 O and 2 H data. Consistent results were obtained with respect to three evaluated elements of lake water balance (I TOT /E ratios, MTT values and the G index) reported in Table 2.These quantities, derived independently from 2 H-based and 18 O-based isotope mass balances are highly correlated (R 2 = 0.9681, 0.9972 and 0.9745 for I TOT /E ratios, MTT values and the G index, respectively).The total inflow-to-evaporation ratios derived from the 18 O-based balance turned out to be ca.10.8 % higher on average than those derived from the 2 H-based balance.For the G index this difference is approximately 2.3 %.The MTT values were ca.12.5 % higher for the 2 H-based balance.
Small but significant differences in the values of the evaluated quantities (I TOT /E ratios, MTT, G index), derived independently from 18 O-and 2 H-based isotope mass balance calculations, stem most probably from their different sensitivity to small changes of the measurable parameters (air and lake water temperature, relative humidity, isotopic composition of local precipitation, isotopic composition of lake water) rooted in the different role of equilibrium and kinetic fractionation during the evaporation process.While for 18 O the ratio of equilibrium to kinetic isotope fractionation is in the order of one, for 2 H it is ten times higher.Since the isotope mass balance method relies on isotope enrichment of lake water along the local evaporation line, controlled mostly by kinetic fractionation, the 18 O-based balance calculations are generally considered more reliable (e.g.Rozanski et al., 2001).

Uncertainty assessment
The above methodology for quantifying elements of water balance in the lakes studied introduces some uncertainties linked to the assumptions made and the uncertainties associated with the parameters used in the evaluation process.Sensitivity tests were performed to derive the range of uncertainties associated with the quantities being evaluated, such as mean turnover time, total inflow-to-evaporation ratio and the G index.The sensitivity analysis was focussing on Eq. ( 7).All variables present in this equation were considered in this process.The results for 18 O-based calculations are summarized in Table 3.
The uncertainty with regard to lake water temperature was probed assuming the temperature change of ±0.87 • C (see Sect. 4.1).The uncertainty of lake water temperature leads to uncertainty as regards the evaporation flux, which in turn influences the I TOT /E , MTT and G values derived for each lake.Also equilibrium isotope enrichment is a function of temperature.The mean turnover time increases by ca. 15 % when the temperature of the lake is reduced by 0.87 • C, and decreases by approximately 12 % when the temperature increases by the same amount.The G index reveals lower sensitivity (2.9 and 3.6 %, respectively).The smallest changes were obtained for I TOT /E ratios (0.9 and 0.7 %, respectively).
The changes of relative humidity of the atmosphere normalized to the temperature of the lake surface have an impact on the evaporation flux, control the I TOT /E ratios through Eq. ( 7) and determine the actual value of kinetic isotope enrichment ε.It was assumed in the calculations that normalized relative humidity changes by ±2 %.As seen in Table 3, the resulting changes of the derived quantities are moderate, the mean turnover time being the most sensitive parameter.
It is apparent from Table 3 that among isotope parameters occurring in Eq. ( 7), the isotopic composition of lake water  (δ LS ) and the isotopic composition of the total inflow (δ IT ), are the two most important variables in the isotope mass balance calculations.An increase of δ 18 O LS by 0.5 ‰, which may account for possible departures from the isotopic steadystate of the investigated lakes (see Sect. 4.4), leads to a decrease of calculated I TOT /E ratios on average by 15.8 %, an increase of MTT values by 19.6 % and a decrease of G index values by 3.3 %.An increase of δ 18 O IT by 0.5 ‰, which may result from the contribution of an evaporated lake water originating from an upstream lake to groundwater input, leads to a substantial increase of I TOT /E ratios (20.2 % on average), a comparable decrease of MTT values (14.8 % on average) and a moderate increase of the G index (4.4% on average).Variation of the 18 O isotopic composition of atmospheric water vapour by ±1.0 ‰ introduces changes in the calculated elements of the water balance of the studied lakes in the order of several per cent (Table 3).
Figure 9 shows the percentage changes of I TOT /E ratios calculated for all studied lakes using Eq. ( 7), in response to the increase of δ LS or δ IT by 0.5 ‰.The sensitivity of the calculated I TOT /E ratios to the given increase of δ LS or δ IT raises sharply with increasing values of this parameter.This is particularly true for through-flow systems characterized by high I TOT /E ratios.Therefore, when isotope studies aimed at quantifying the water balance of such systems are planned, it is important to characterize these two isotope quantities as accurately as reasonably possible.

Conclusions
The Rokua esker, with its numerous lakes located across a relatively small area, provided a unique opportunity to explore the possibilities offered by environmental isotope techniques in quantifying the water balances of lakes and their dependency on groundwater in a sub-polar climatic setting.The quantification of groundwater seepages to lakes using conventional methods is notoriously difficult and associated with considerable uncertainty.The presented study demonstrates the power of isotope mass balance approach for resolving this issue.It appears that a stable isotope analysis of lake water samples, collected at the right time and supplemented by appropriate field observations, may lead to quantitative assessment of the water balance of a large number of lakes located in a similar climatic setting.
The presented study has demonstrated that consistent results are obtained when the isotope mass balance calculations are run independently for oxygen-18 and deuterium.This strengthens the position of heavy stable isotopes of water as a unique tool for quantifying elements of the water bal-E.Isokangas et al.: Quantifying groundwater dependence of a sub-polar lake cluster ance of lakes, particularly for groundwater-dominated systems.Solving three equations simultaneously (one water balance equation plus two isotope balance equations) may help to quantify key balance-related parameters such as evaporation and groundwater inflow and outflow rates for the studied lake system, which are difficult to quantify using conventional methods.
The specific behaviour of lakes located in sub-polar regions, with their seasonal ice cover extending over several months, offers another opportunity for quantifying groundwater seepage during ice-cover periods.As shown in this study, observations of seasonal changes in the stable isotopic composition of lake water, in particular during the ice-cover period, combined with the survey of isotopic composition of groundwater in the vicinity of the lakes studied, allows the quantification of groundwater fluxes to this lake during winter.If such an approach is combined with the isotope mass balance calculations performed for the ice-free summer season, important information about the seasonal variability of groundwater seepage to lakes located in sub-polar and polar regions can be obtained.
The G index characterizing the groundwater dependency of a lake proposed in this study, and defined as a percentage contribution of groundwater inflow to the total inflow of water to the given lake, appears to be a straightforward, quantitative measure of this dependency.The studied Rokua lakes appear to be strongly dependent on groundwater; more than 40 % of water received by these lakes comes as groundwater inflow.The quantitative evaluation of groundwater dependency of lakes via the G index proposed in this study may assist lake restoration policies in areas where groundwater is a source of nutrients to the studied lakes.

Figure 2 .
Figure2.The long-term (1959The long-term ( -2013) )  monthly mean values of surface air temperature and the amount of precipitation recorded at the station located 10 km south-west of the study site (Finnish Meteorological Institute, 2014).

Figure 4 .
Figure 4. δ 2 H-δ 18 O relationship for different appearances of surface water (lakes, streams) and groundwater in the study area, investigated within the scope of this study.The Rokua evaporation line (local evaporation line -LEL) was defined as the best fit line of the data points representing lakes sampled during the July-August 2013 campaign.

Figure 6 .
Figure 6.Seasonal variations of δ 18 O in lake Ahveroinen 1, observed at 1 and 4 m depth.Maximum depth of the lake is 4.8 m.

Figure 7 .
Figure 7.The ratio of total inflow to evaporation (I TOT /E) as a function of the measured 18 O isotope enrichment ( δ 18 O) for Rokua lakes surveyed during the July-August 2013 sampling campaign.

Figure 8 .
Figure 8.The G index quantifying the groundwater dependency of lakes in the Rokua study area.The index is defined as the percentage contribution of groundwater inflow to the total inflow of water to the given lake.Shown are the mean G values obtained from independent isotope mass balance calculations based on 2 H and 18 O data.
The (+) and (−) signs signify an increase or reduction, respectively, of the derived quantity by the reported percentage.

Figure 9 .
Figure 9. Changes of the total inflow-to-evaporation ratio (in %) based on 18 O isotope mass balance as a function of I TOT /E value, in response to parameter change, calculated for the studied lakes on the Rokua esker.Two cases are considered: (a) an increase of the measured δ 18 O of lake water by 0.5 ‰, and (b) an increase of δ 18 O of the total inflow to the given lake by 0.5 ‰.See text for details.

Table 1 .
Mean isotopic composition of selected lakes, streams and groundwater, sampled four times per year during the period 2010-2012.

Table 2 .
The results of isotope mass balance calculations for 67 lakes of Rokua esker sampled during the July-August 2013 field survey.

Table 3 .
Sensitivity of selected elements of the 18 O-based water balance of the studied lakes to changes of the parameters involved.