Temporal and lake-specific variations in oxygen and hydrogen stable isotopes in a boreal lake-chain during two hydrologically differing years

ABSTRACT

Hydrogen and oxygen stable isotope mass ratios (δ²H and δ¹⁸O) are conservative tracers that reflect hydrological conditions of watersheds and lakes, information that can be used for water budgets, such as estimating nutrient loading from the catchments and direct and indirect effects of climate change on waterbodies. We measured δ²H and δ¹⁸O values and estimated water mass balances over 2 hydrologically differing years in a chain of 10 boreal lakes. During a warm and dry summer, evaporative loss enriched the isotope values at lake surfaces. Thus, in lakes that were thermally stratified during summer, the isotope values differed distinctively between epilimnion and hypolimnion. During the following rainy summer, the within-lake variation of isotope values was much lower with no clear trend, presumably a result of a more unstable water column and faster water throughflow in the lake chain. After autumnal mixing the lake δ¹⁸O values were strongly correlated with lake area and proportion of inland waters of the catchment area. The interannual within-lake difference was then ∼1‰, indicating differences in evaporation. Isotope mass balance (IMB) calculations and those based on lake temperature and hydrological measurements in the nearby region yielded similar results for the rainy ice-free period. Presumably, prerequisites of steady-state conditions for IMB were better met in the rainy year than during the previous drier year. For some lakes the results suggested remarkable volumes of groundwater inflow, whose role in boreal lake water budgets deserves more attention in future studies.

KEYWORDS:

• boreal lake-chain
• dry/wet year
• hydrogen-2
• oxygen-18
• precipitation; stable isotopes

Introduction

Stable isotope mass ratios of hydrogen (δ²H) and oxygen (δ¹⁸O) and the knowledge of kinetic fractionation against heavier molecules in water reactions are widely applied in hydrological research to estimate water origin, evaporation, and flow through a landscape, starting from precipitation and ending through catchments to rivers and lakes (e.g., Gibson 2002, Gibson and Edwards 2002, Jonsson et al. 2009, Haig et al. 2021). Global variation in precipitation and its original isotope mass ratios as well as seasonal variation in the ratio between evaporation to precipitation are strongly reflected in the δ²H and δ¹⁸O values in soils and aquatic systems (Kortelainen and Karhu 2004, Gibson et al. 2008, Jonsson et al. 2009). Stable isotope values are also used to model groundwater flow to lakes and rivers (Stets et al. 2010, Isokangas et al. 2015), but few lake studies have covered time series of δ²H and δ¹⁸O values over different seasons and depth profiles, including both epilimnion and hypolimnion. Among those, Halder et al. (2013) applied an isotope approach for studying transport and mixing of River Rhone water in Lake Geneva in Switzerland and France; Perini et al. (2009) studied mixing, stratification, and water origin in 6 lakes in the Italian Alps; and Braig et al. (2010) investigated these phenomena in 3 lakes in southern Germany.

Knowledge on progressive enrichment of heavier isotopes during evaporation in waterbodies has been used to estimate regional water balances, especially in remote areas where direct hydrological measurements are infrequent and/or sparse (e.g., Gibson and Edwards 2002, Bennett et al. 2008, Fellman et al. 2011, Shi et al. 2017, Yapiyev et al. 2020). The isotope mass balance (IMB) methods are based on the model originally described by Craig and Gordon (1965) and further revised and applied by, for example, Gat et al. (1994), Horita and Wesolowski (1994), Gibson (2002), Bennett et al. (2008), Horita et al. (2008), Gibson and Reid (2014), and Gibson et al. (2016). In addition to measurements of δ²H and δ¹⁸O in the waterbodies, respective values for atmospheric moisture are needed, but these are seldom directly measured; usually, the values for rainwater collected from the region or derived from global data sources are used (e.g., Bennett et al. 2008, Yapiyev et al. 2020). The complexity of calculations and the sensitivity of results with δ²H and δ¹⁸O values of ambient air moisture have proved challenging for a wider application of the IMB method. Skrzypek et al. (2015) tested the effects of evaporative loss on δ²H and δ¹⁸O enrichment in pan measurements in arid conditions and successfully applied the equations of the previously noted papers in the IMB models, combined in Hydrocalculator software, making the applications more user friendly for people not familiar with isotope chemistry.

We studied δ²H and δ¹⁸O values of water in a lake chain consisting of 10 lakes in southern Finland. The study lakes differ in their topographic, morphometric, and landscape properties. All study sites belong to the same large water system, which means that water from all the lakes flows through the lowermost and largest lake in this hydrological pathway. Therefore, we were able to follow how the stable isotope signatures were changing along the pathway from the upstream to downstream lakes. We used both hydrological and limnological measurements in the study lakes and in the nearby region as well as stable isotope ratios for water balance calculations in the lake chain. This study is a part of a project investigating loading, retention, and metabolism of carbon, nitrogen, and phosphorus in the upper reaches of a large boreal water basin, River Kokemäenjoki, running to the Baltic Sea (Valpola and Ojala 2006). Because the 2-year study period (2003 and 2004) had extraordinary hydrological conditions, with the lowest groundwater table during the first year and the highest summer flood during the second year since the early 1970s, the hydrological properties of the lakes were especially emphasized. We focused on 3 major questions: (1) how do seasonal and interannual hydrological conditions affect water δ²H and δ¹⁸O mass ratio values, (2) how do morphometric properties of the lakes and their position in the landscape affect these ratio values, and (3) do water balance estimates based on physical measurements in the lakes and nearby regions yield similar results to those based on stable water isotopes?

Background information

Study area

The study lakes form a hydrological network within a catchment area of 1473 km² in southern Finland. The area is in the boreal zone, situated in the uppermost reaches of the Kokemäenjoki watercourse, the fifth largest drainage system in Finland, and flows to the Baltic Sea (Fig. 1, Supplemental Table S1). The landscape lies on granodiorite and gneiss bedrock, and glacial till is the dominant soil type with podzol soil on top (Valpola and Ojala 2006). Locally glaciofluvial soils (granules, sands, and silts) influence hydrology by groundwater formation (Valpola and Ojala 2006, Arvola et al. 2015). The catchment areas of the lakes are mainly covered by forests (40–71%), which are largely managed by forestry operations (logging, ditching). About 80–90% of peatlands in southern Finland have been ditched during the past 60 years and are included in the forests (Supplemental Table S1; data for catchment area and characteristics obtained by VALUE tool, http://paikkatieto.ymparisto.fi/value/). Agricultural land and urban environments (villages and industrial area) cover 7–35% and 2–21%, respectively, of the catchments of the study lakes while the open wetlands only cover a small proportion (0.1–4.2%).

The study lakes are typical headwater lakes in southern Finland by means of lake size (area 0.05–31.5 km²), morphology (maximum depth 3–73 m), water quality, and human influence. The lakes are connected by rivers and streams within ∼30 km from the inflow to the uppermost lake down to the outflow of the lowermost lake. The elevation difference between the highest and lowest site is 9.9 m at the mean water table. The lowest lake in the lake chain, Lake Iso-Roine, is also interconnected via River Vihavuosi with another branch of several upstream lakes not included in the study. The largest lakes (Ormajärvi, Suolijärvi, Pyhäjärvi, and Iso-Roine) are dimictic, whereas 3 shallow lakes (Lehee, Kyynärö, and Teuronjärvi) are polymictic. Two small headwater lakes, Lampellonjärvi and Lovonjärvi, are meromictic with permanent anoxic hypolimnion (data not shown). No water column data are available for Lake Takanen.

The groundwater sampling site, a local municipal raw water pumping station, is situated by the glaciofluvial sand ridge next to Lake Ormajärvi (Fig. 1), thus representing the local groundwater characteristics.

Figure 1. Catchment of the study area (black line) with subcatchments (red lines) and study lakes (blue) with connecting rivers and their sampling sites (dots). River sites: a = Ormajärvi, outflow; b = Suolijärvi, inflow; c = Suolijärvi, outflow; d = Lehee, inflow; e = Lehee, outflow; f = Pyhäjärvi, inflow; g = Iso-Roine, inflow southern; h = Iso-Roine, inflow northern; i = Iso-Roine, outflow; j = river Alvettulanjoki. WS = weather station, GW = groundwater sampling site, Q = discharge measurement sites. Arrows indicate flow direction.

Climate in the region

Our study area belongs to the northern temperate climate zone where temperature has a distinct seasonal pattern (Jylhä et al. 2014). The long-term (1981–2010) mean annual air temperature from April to October is 4.2 °C and mean monthly temperature typically >0 °C (Fig. 2). Annual precipitation during the same period varied between 500 and 715 mm (average 645 mm). Monthly mean of relative humidity varied between 68% and 93% (measurements from 2002–2010), with lowest values in April and May (all measurements from Lammi Pappila weather station of Finnish Meteorological Institute (FMI, https://www.ilmatieteenlaitos.fi/).

Figure 2. Monthly mean air temperature (T, lines) and precipitation sum (p, columns) during the study period at Lammi Biological Station (called Pappila Meteorological Station in FMI database).

In spring, rapid snowmelt is followed by high inflow to lakes just before ice-out. Larger lakes in the study area are ice-free from April–May to middle November–late January, and small lakes from April to early November–late December, thus ∼1 month shorter than the large lakes. In 2003 and 2004, the observed ice-out and ice-on dates in Lake Pääjärvi, a nearby lake with ice records, were 8 May and 21 December, and 27 April and 13 December, respectively.

Material and methods

Sampling

Five lakes in the chain (Ormajärvi, Suolijärvi, Lehee, Pyhäjärvi, Iso-Roine; Fig. 1, Supplemental Table S1) and river sites between those were sampled 16–17 times in 2003 and 17–18 times in 2004. The lake and river sites were sampled twice a month during the ice-free period (May–Oct), 2–3 times during the late fall before ice cover (Nov–Dec), and 1–4 times during the ice-cover period (Jan–early Apr). In addition, 5 small upstream lakes (Lampellonjärvi, Lovonjärvi, Kyynärö, Teuronjärvi, and Takanen) were sampled once in 2003 (late Oct) and 4 times in 2004 (late Mar, early May, late Aug, late Oct). Lake samples were taken from the deepest point of each lake while river samples were taken at outlets (50–100 m downstream), just before the inlets of next lakes, thus representing incoming and outflowing water quality (sites marked a–j in Fig. 1). In short rivers, samples were collected from one site only. At each site and sampling, dissolved oxygen (DO) concentration and water temperature were measured manually with a YSI ProODO optical dissolved oxygen probe (ysi.com, Yellow Springs, OH, USA) at 1 m intervals (lakes) and in rivers from the surface.

Lakes were sampled 1 m below the surface (i.e., epilimnion) and 1 m above the bottom (i.e., hypolimnion), whereas rivers were sampled only from the uppermost 0.3 m depth to avoid disturbance with sediment. All water samples were taken by a 2 L tube sampler (Limnos) and collected into 20 mL plastic (HD-PE) or glass liquid scintillation vials, which were filled, closed tightly, and stored in a dark and cool (4 °C) environment until analysed within 3–5 days. In addition, we applied long-term δ¹⁸O and δ²H results of precipitation collected at local meteorological stations (Pappila during 1998–2003 and Evo station, ∼20 km northward, during 2004–2010). The samples were collected monthly into 1 L HD-PE bottles with a funnel and paraffin oil layer to prevent evaporation. Groundwater was sampled 3 times (May, Jul, and Sep 2004) at the local municipal raw water pumping station.

Isotope analyses

Water samples of 0.5 mL were equilibrated with CO₂ at 25 °C for oxygen isotope analysis and 0.3 mL with H₂, Pt as catalyst, at 25 °C for hydrogen isotope analysis. Equilibrated gas was introduced online via an interface (Gas Bench II) to a mass spectrometer (ThermoFinnigan Delta^Plus XL, Thermo Scientific, Bremen, Germany) to measure isotope ratios. Analyses were conducted at Laboratory of Chronology, University of Helsinki, Finland.

The δ¹⁸O and δ²H values are expressed as δ notation in per mil relative to Vienna Standard Mean Ocean Water (VSMOW) standard: (1) ${\delta }^{18}\mathrm{O}(\mathrm{or}{\delta }^2\mathrm{H})\text{‰}=\left({\displaystyle \frac{\mathrm{R}(\mathrm{sa})}{\mathrm{R}(\mathrm{st})}-1}\right)\times 1000,$(1) where R is a ratio of ¹⁸O/¹⁶O or ²H/¹H for a sample (sa) and a standard (st), respectively.

The laboratory reference samples were normalized on the VSMOW–SLAP scale (Coplen 1988). Errors for the δ¹⁸O and δ²H values were ±0.15‰ and ±1‰, respectively. Because the number and accuracy of measurements for δ²H values were lower, the main focus of our results is on changes in δ¹⁸O values in the lake chain.

Water balance

We estimated water balance in the study region using (1) data for hydrological variables (precipitation, air and water temperature, relative humidity, and discharge) measured in each lake and at the nearest monitoring sites of the study area, and (2) the IMB method using measured δ¹⁸O and δ²H values in the study lakes, their inflows, and precipitation.

Hydrological measurements

Data for daily precipitation were obtained from the open database of FMI (https://en.ilmatieteenlaitos.fi/). We use here data from 4 weather stations located ∼5–40 km from the study lakes. Data for catchment runoff (4 stations), discharge (4 stations), and evaporation (pan measurements at 1 station) measured close to the study region during the study years were taken from Hertta open database of Finnish Environment Institute (https://www.syke.fi/avoindata; Table 1). Discharge from each subcatchment to the study lakes was estimated using daily measurements for 2003 and 2004 from rivers Mustajoki (for the 5 small upstream lakes), Teuronjoki (for Ormajärvi, Suolijärvi, Lehee, and Pyhäjärvi), and Vihavuosi (for Iso-Roine), all situated <30 km from the study sites and having land use characteristics and approximate size similar to those of the study lakes. We assumed that the discharge (Q_a; m³ s⁻¹) from a subcatchment (CA_a) was relative to its surface area compared to the measured values (Q_b; m³ s⁻¹) from the nearby respective subcatchment (CA_b; i.e., Mustajoki, Teuronjoki, and Vihavuosi discussed earlier): (2) ${\mathrm{Q}}_{\mathrm{a}}={\displaystyle \frac{\mathrm{C}{\mathrm{A}}_{\mathrm{a}}}{\mathrm{C}{\mathrm{A}}_{\mathrm{b}}}\times {\mathrm{Q}}_{\mathrm{b}}.}$(2)

Table 1. Precipitation (average and SD of 4 meteorological stations: Pirttikoski, Pappila, Iso-Evo, Palomaa) in 2003 and 2004; data from open database of Finnish Meteorological Institute (FMI). Runoff from the catchment (average and SD of 4 small brooks and rivers: Valkea-Kotinen, Löytynoja, Mustajoki, Teuronjoki), discharge (average and SD) in 3 rivers (catchment area in brackets) and evaporation (Jun–Sep) measured at Jokioinen FMI station (pan measurements) during the study years. Hydrological data are from https://www.syke.fi/avoindata, an open database of Finnish Environment Institute. All stations are located at 30–90 km from the lake chain.

	2003		2004
	average	SD	average	SD
Precipitation sum; mm
Whole year	596.4	19.4	772.6	24.6
May–Oct	393.1	20.6	553.2	27.9
Jun–Sep	196.1	10.6	460.2	20.4
Runoff; L s^−1 km^−2
Whole year	5.47	1.64	12.10	3.03
May–Oct	6.21	1.64	15.01	3.86
Jun–Sep	2.98	1.56	17.6	4.50
Discharge; m^3 s^−1
Mustajoki (84 km^2)	0.34	0.49	0.78	1.06
Teuronjoki (244 km^2)	1.10	1.19	2.97	2.17
Vihavuosi (863 km^2)	3.07	1.32	8.57	2.32
Evaporation Jokioinen sum; mm
Jun–Sep	386.7		314.2

The calculated Q_a accounts for evaporation for subcatchment CA_a in similar proportion to that in the measured discharge (Q_b) in the river collecting waters from subcatchment CA_b because of their similar characteristics (fields, forests and upstream lakes). Total water inflow (I) to each lake during the ice-free period (May–Oct) was estimated by adding the volume of precipitation (P) on the lake surface during the same period: (3) $\mathrm{I}=\mathrm{Q}+\mathrm{P}.$(3)

Evaporation (E; mm) from lake surface was calculated with the bulk aerodynamic method using equation by Järvinen and Huttula (1982): (4) $\mathrm{E}=(0.236+0.068V_2)(e_0-e_2),$(4) where V₂ is the wind velocity 2 m above ground (m s⁻¹), e₀ is the saturation vapor pressure (mb) at lake surface temperature (1 m) and air humidity, and e₂ is the prevailing vapor pressure (mb) at 2 m height. Data for wind velocity, air temperature, and humidity were taken from hourly measurements at Jokioinen meteorological station (∼70–80 km southwest of the study area). Because lake water temperature was not intensively measured in the 5 small upstream lakes, we assumed evaporation in these lakes was of similar magnitude and thus used the obtained mean May–October E values for those estimates.

Isotope mass balance (IMB)

The IMB balance model assumes progressive enrichment of heavier isotopes due to evaporation in the waterbodies downstream of the drainage area, the steady state of the lakes, and the evaporation (E) in a lake equaling or exceeding the inflow (I) from the upstream lakes and rivers (Gibson and Reid 2014, Skrzypek et al. 2015). Thus, we used data collected in both years in October from all 10 lakes and mean values for the ice-free periods (May–Oct), which included δ¹⁸O measured 11 times from 5 lakes in 2003 and 2004 (Ormajärvi, Suolijärvi, Lehee, Pyhäjärvi, Iso-Roine) and 3 times from the small headwater lakes in 2004. We used the isotope values for samples from 1 m depth for IMB calculations, although most lakes were stratified in June–September (Supplemental Fig. S1). We assumed that lakes with the position 1–2 in the landscape (Lampellonjärvi, Lovonjärvi, Kyynärö; Fig. 1) received mainly groundwater. For the inflow to Ormajärvi, Teuronjärvi, and Takanen (position 3 in the lake chain), we used mean isotope values of those with the position 1–2 in the landscape. For Suolijärvi inflow we used the catchment area weighted isotope values of Ormajärvi and Teuronjärvi. For Lehee inflow we used isotope values for the inflow from the lake’s main inflow, Suolijärvi. For Pyhäjärvi we used inflow values from Lehee, and for Iso-Roine we used catchment area weighted values from Pyhäjärvi and Vihavuosi outflows. We used mean δ¹⁸O values of rainwater collected at 2 nearby stations (Fig. 1) during May–October in 2003 and 2004 for the values of atmospheric moisture. For a comparison we used δ²H and δ¹⁸O values of precipitation collected from Espoo (close to the Baltic Sea, ∼115 km south of the study area) for the respective periods from the database of Global Network of Isotopes in Precipitation (GNIP, International Atomic Energy Agency [IAEA]; https://nucleus.iaea.org/wiser/index.aspx). However, because the δ²H and δ¹⁸O values of precipitation in Espoo were significantly more enriched than those measured in Lammi (1.6‰ and 12.1‰, paired t-tests p < 0.01 for both, df = 19), we used the values from Lammi for the IMB estimates. For all calculations we used average air temperature and relative humidity (May–Oct) measured hourly at the nearby Lammi Pappila weather station. Evaporation to inflow ratio (E/I) was calculated with Hydrocalculator software (Skrzypek et al. 2015): (5) ${\displaystyle \frac{E}{\mathrm{I}}=\left[{\displaystyle \frac{({\delta }_L-{\delta }_P)}{({\delta }^{\ast }-{\delta }_L)\times m}}\right],}$(5) where δ_L is the δ¹⁸O value of the studied lake and δ_P is the respective value of the lake or inflow river above, δ* is the threshold of limiting isotopic composition, and m is the enrichment slope. The latter 2 parameters were calculated by Hydrocalculator (Supplemental Tables S2, S3). Local evaporation line (LEL) was calculated from the regression slope of δ²H and δ¹⁸O values of the intensively studied 5 lakes (1 m) during both ice-free seasons. We also tried step-by-step calculation of E/I using biweekly obtained δ¹⁸O values from the 5 lakes, but these yielded unrealistic estimates, presumably because of the uncertainty of ambient air moisture δ¹⁸O values based of rain water during the same period (data not shown).

For the warm and dry summer period from June to August 2003 we also applied the non-steady state model: (6) $f=1-{\left[{\displaystyle \frac{({\delta }_{\mathrm{La}}-{\delta }^{\ast })}{({\delta }_{\mathrm{Lj}}-{\delta }^{\ast })}}\right]}^{{\displaystyle {}_{\frac{1}{m}}}},$(6) where f is the evaporative loss fraction of the lake volume, δ_La is the δ¹⁸O value on 18 August, and δ_Lj the value on 16 June, and the parameters δ* and m are as in equation 4.

Other limnological parameters

Schmidt stability (S) of the water column (J m⁻²) in the lakes (Schmidt 1928, Idso 1973) was calculated with the equation presented by Korhonen (2002): (7) $\mathrm{S}={\displaystyle \frac{\mathrm{g}}{{\mathrm{A}}_0}\underset{0}{\overset{\mathrm{h}}{\int }}\mathrm{A}(z_{\mathrm{g}}-z)({\rho }_0-\rho ){\mathrm{d}}_{\mathrm{z}},}$(7) integrated from lake surface to bottom (0 to h), where A₀ is the lake surface area; A is lake area at the depth z; g is acceleration due to gravity; z_g is the depth of geometric gravity point of the lake; ρ₀ and ρ are water densities at lake surface and at the depth z, based on temperature profiles; and d_z is total depth of the water column.

The hydraulic residence time of each lake was estimated by dividing lake volume by annual surface water inflow (cf. equation 3). For stratified lakes, residence time was also estimated from epilimnetic volume, based on the vertical distribution pattern of water temperature and dissolved oxygen showing the depth layer of the greatest difference of these within 1 m distance.

Results

Hydrological conditions

The study period until late autumn 2003 was exceptionally dry from a drought that began during summer 2002 (precipitation ∼80% of the long-term average). Because of the low precipitation, the groundwater table was exceptionally low during 2003 (data not shown; see https://www.syke.fi/avoindata). During the first study year in 2003, summer (Jun–Aug) mean air temperature at Lammi Biological Station was 15.6 °C and in 2004 was 14.4 °C; maximum air temperature was 30.7 and 27.6 °C, respectively. Correspondingly, in July 2003 mean air temperature was 19.6 °C and relative humidity (Rh) 77.2%, and a year later the values were 15.7 °C and 82.4%. In agreement with the Rh values, June–July precipitation was 3 times higher in 2004 than in 2003 (i.e., 326 and 107 mm; Fig. 2). In 2003, total precipitation in the study area was 596 mm but during the following year was ∼30% higher (average of 4 FMI stations; Table 1). Respectively, river discharge values were twice as high in 2004 than in 2003. For example, in River Mustajoki, a nearby small river with an online hydrological station, annual average runoff in 2004 was the second highest since the onset of the 1971 measurements. A similar increase in discharge took place in River Vihavuosi, the largest inflow to Lake Iso-Roine (Table 1, Fig. 1).

The differences in the weather and lake morphology affected the stability of water columns of the study lakes (Supplemental Fig. S1, Supplemental Table S4). In 2003, maximum water temperature at 1 m depth varied between 22.4 and 24.5 °C, whereas in 2004 the range was between 19.6 and 20.5 °C. In Lake Lehee, water column stability was low during both years. In other intensive study lakes, water column stability from early June to the end of August was 45–76% lower in 2004 than in 2003.

Seasonal and interannual variation in δ²H and δ¹⁸O

Annual long-term mean stable isotope mass ratios in precipitation, δ²H = −89.8‰ and δ¹⁸O = −12.6‰, were close to those measured in groundwater in 2004 (−92.2‰ and −12.5‰, respectively). During the 2 successive years, the δ²H and δ¹⁸O isotope values of precipitation followed a similar seasonal pattern, although in 2003 values were slightly lower than in 2004 (Fig. 3, Supplemental Table S5). Seasonally, lower values were recorded in winter than in summer, as seen in the long-term (1981–2010) data. Local meteoric water line (LMWL y = 7.575x + 5.350; c.f. Rozanski et al. 1993, Putman et al. 2019) is based on the linear relationship between δ²H and δ¹⁸O isotope values in this long-term dataset, together with the values for precipitation collected from Espoo during the same period (Fig. 4, Supplemental Fig. S2). The slope of the LEL was 5.45 based on isotope mass ratios of δ²H and δ¹⁸O in the 5 intensively studied lakes during both years (Supplemental Fig. S3).

Figure 3. Oxygen isotope ratio (δ¹⁸O‰) in precipitation (average with range) collected monthly at Pappila (1998–2003) and Evo (2004–2010) FMI meteorological stations. The distance between stations is 19 km. The values for the study years are marked as open triangles (2003) and black dots (2004).

Figure 4. Local meteoric water line (straight line) is based on δ²H and δ¹⁸O measurements of precipitation collected monthly at Pappila (1998–2003) and Evo (2004–2010) meteorological stations as well as from Espoo, Finland (GNIP database; Supplemental Fig. S2). Local evaporation line (hatched line) is based on epilimnetic values of intensive study lakes during both study years. Mean δ²H and δ¹⁸O values for precipitation (Mean prec) during 1998–2010 and groundwater (GW) measured at one site in 2004 (n = 3) are also shown (overlapping dot and circle). VSMOW = Vienna Standard Mean Ocean Water.

The δ¹⁸O and δ²H values in the uppermost headwater lakes were closer to those in precipitation and groundwater than values from lakes lower in the chain (Supplemental Table S5). Except for Lake Suolijärvi, the δ¹⁸O and δ²H values gradually enriched during both study years along the lake chain. The most enriched δ¹⁸O and δ²H values were measured at the Alvettulanjoki River sampling site (Fig. 1), which also received water from the Hauhonselkä catchment area (1585 km²). At this site, δ¹⁸O was 2.4–3.1‰ and δ²H 13.8–19.9‰ more enriched than those in the epilimnion of the lakes with the position 1–2 in the chain.

In the lakes, seasonal variation in the isotope values was greater in 2003 than in 2004 (Fig. 5, Supplemental Table S5). In 2003, the epilimnetic δ¹⁸O values in the 5 intensively studied lakes started to increase during the summer season and resulted in a clear deviation (0.5–1‰) relative to those in the hypolimnion. The only exception was the polymictic and shallow Lake Lehee, where the whole water column was homogeneous. The enrichment of the δ¹⁸O values in the epilimnion was linear from early June to late August. In summer 2004, the epilimnetic isotope values did not increase in the study lakes except in Lake Suolijärvi, where the values slightly increased, separating them from the hypolimnion. At the time of complete autumnal mixing, the difference in isotope values between the epilimnion and hypolimnion vanished in all lakes.

Figure 5. Time series of the δ¹⁸O values (‰) at 1 m depth and 1 m above bottom of 5 study lakes from March 2003 to November 2004. Circled dots denote samples taken under ice.

In winter, beneath the ice cover, the δ¹⁸O isotope values in Ormajärvi 0.3–0.5‰ were lower close the bottom than at 1 m depth (both years), whereas in Pyhäjärvi more enriched values were found close the bottom during winter 2004. In the other intensive study lakes, no clear difference was found between the 2 layers during winter. The only exception was a low value (12.4‰) measured on 4 April 2003 in Lake Suolijärvi when snowmelt and rainwater accumulated on the ice surface, presumably contaminating the water sampled from 1 m depth.

The epilimnetic δ¹⁸O values measured in October during both study years were strongly correlated with lake area, proportion of inland waters in the catchment area and total cumulative catchment area, and lake position in the landscape, all intercorrelated variables (Fig. 6, Supplemental Table S6). During both study years, the δ¹⁸O enrichment in the lake chain followed the same pattern, but during the wet year (2004) the values were ∼1‰ lower, and the intercept of the regression lines of the 2 years is significantly different (ANCOVA, t = 4.044, p < 0.001, df = 19).

Figure 6. δ¹⁸O values (‰) of epilimnetic waters in all 10 study lakes in October 2003 and 2004, related to (a) lake area and (b) proportion of inland waters of the catchment, including the lake itself. The intercepts of the regression lines differ significantly (ANCOVA, p < 0.05). The circled values in (a) refer to Lake Kyynärö, which is presumably more influenced by groundwaters than the other lakes.

Water balance

The interannual differences in precipitation and temperature were reflected in water balances of the study lakes (Table 2). During the 2004 rainy year, the water inflow to the lakes was >2 times higher than during the previous year. Evaporation (E) from lake surfaces during May–October 2003, calculated with bulk aerodynamic method, varied between 2.6 and 3.0 mm d⁻¹ (sum 505–544 mm), ∼30% higher than precipitation during the same period (Table 1). The average lake evaporation was 24% higher in 2003 than during the open-water season of the following year (2.2–2.4 mm d⁻¹) when precipitation exceeded evaporation. A similar interannual difference was observed at Jokioinen, where the Finnish Meteorological Institute conducts pan measurements (Table 1; data only for Jun–Sep). The estimated theoretical residence time of the study lakes was highly variable, ranging from ∼0.08 yr in Lake Lehee in 2003 to 10 yr in Lake Ormajärvi (Supplemental Table S6). When considering only epilimnetic volume in thermally stratified lakes, the residence times were much shorter (1.9–3.6 yr) than the theoretical estimates. During the 2004 rainy year, the residence times were about one-third lower than during the previous year and had similar between-lake variation.

Table 2. Lake area (LA), catchment area (CA), position in the landscape (Po), and inflow (I) estimated on catchment area-related inflow measurements of respective size catchments (equations 2–3), evaporation from lake surface (E) based on bulk aerodynamic calculations (equation 4) from May to October in 2003 and 2004. Whole lake evaporation estimates for the 5 small upstream lakes* used average E estimates (2.8 and 2.3 mm d⁻¹ for 2003 and 2004, respectively) of the 5 intensively studied lakes. mm = 0.001 m³ m⁻², nd = no data.

Lake	LA	CA	Po	I (May–Oct) 10^3 m^3		E (May–Oct) mean mm d^−1		E (May–Oct) sum mm		E (May–Oct) 10^3 m^3
	km^2	km^2		2003	2004	2003	2004	2003	2004	2003	2004
Lampellonjärvi	0.06	1.3	1	109	262	nd	nd	nd	nd	31*	25*
Lovonjärvi	0.05	6.7	2	476	1253	nd	nd	nd	nd	26*	21*
Kyynäröjärvi	0.25	28.3	2	2022	5305	nd	nd	nd	nd	128*	104*
Teuronjärvi	1.34	39.7	3	3158	7807	nd	nd	nd	nd	685*	557*
Takanen	0.37	16.4	3	1245	3158	nd	nd	nd	nd	189	154
Ormajärvi	6.53	74.8	3	8334	20 655	3.0	2.2	544	413	3605	2643
Suolijärvi	2.05	183	4	16 093	46 305	2.7	2.2	505	405	1018	830
Lehee	1.05	215	7	18 487	53 990	2.6	2.2	526	412	502	425
Pyhäjärvi	9.48	294	8	27 763	76 270	2.8	2.3	509	421	4884	4012
Iso-Roine	31.5	1473	14	108 645	284 326	2.8	2.4	507	440	16 229	13 910

During May–October 2003, the ratio of evaporation to surface waterflow (E/I), estimated with bulk aerodynamic calculations and catchment area ratios, varied between 0.03 and 0.43 in the study lakes, whereas the IMB estimates varied between 0 and 0.16, and even negative estimates were obtained for lakes Suolijärvi and Iso-Roine (Table 3). For the respective period in 2004, the E/I estimates with the different methods were much closer and correlated with each other (Supplemental Fig. S4), varying between 0.01 and 0.16. The IMB estimates based on mean May–October δ¹⁸O values were close or identical to those based only on October values. The E/I estimates adjusted with LEL were ∼1–7% lower than those without the adjustment.

Table 3. Evaporation to inflow ratio estimates for the lakes during open water seasons (May–Oct 2003 and 2004) based on bulk aerodynamic (B) calculations and catchment area (CA) related inflow (Table 2). Respective isotope mass balance results estimated using δ¹⁸O values of the samples taken in October and mean values at 1 m depth during May–October, calculated with Hydrocalculator without and with local evaporation line (LEL) adjustment (see equation 5 and Skrzypek et al. 2015).

	2003					2004
	B&CA based	δ^18O (no LEL adj)	δ^18O (LEL adj)	δ^18O (no LEL adj)	δ^18O (LEL adj)	B&CA based	δ^18O (no LEL adj)	δ^18O (LEL adj)	δ^18O (no LEL adj)	δ^18O (LEL adj)
Lake	May–Oct	Oct	Oct	May–Oct	May–Oct	May–Oct	Oct	Oct	May–Oct	May–Oct
Lampellonjärvi	0.28	0.14	0.10	nd	nd	0.10	0.06	0.04	0.06	0.04
Lovonjärvi	0.05	0.14	0.10	nd	nd	0.02	0.04	0.03	0.04	0.02
Kyynäröjärvi	0.06	0.08	0.06	nd	nd	0.02	0.04	0.03	0.03	0.02
Teuronjärvi	0.22	0.21	0.12	nd	nd	0.07	0.05	0.03	0.04	0.02
Takanen	0.15	0.16	0.09	nd	nd	0.05	0.05	0.03	0.04	0.03
Ormajärvi	0.43	0.15	0.09	0.14	0.08	0.16	0.11	0.06	0.11	0.07
Suolijärvi	0.06	−0.05	−0.03	−0.03	−0.02	0.02	0.01	0.01	0.05	0.03
Lehee	0.03	0.00	0.00	0.04	0.02	0.01	0.00	0.00	0.01	0.00
Pyhäjärvi	0.18	0.01	0.01	0.00	0.00	0.05	0.05	0.03	0.06	0.03
Iso-Roine	0.15	−0.02	−0.01	−0.01	−0.01	0.05	0.01	0.00	0.01	0.00
Average	0.16	0.08	0.05	0.03	0.02	0.05	0.04	0.03	0.04	0.03

During the dry and warm period from 16 June to 18 October, the δ¹⁸O values at 1 m depth in the 5 intensively studied lakes increased linearly (Supplemental Fig. S5). The slopes of the regression lines for this period were similar in the 3 largest lakes, indicating a 0.009‰ daily increase, whereas a greater enrichment occurred in the smaller lakes Lehee and Suolijärvi (0.024‰ and 0.014‰, respectively). The evaporative loss from these lakes was estimated to vary 2–16% (fraction 0.02–0.16) of the lake average volume (Table 4). For the large lakes (Ormajärvi, Pyhäjärvi, and Iso-Roine), the difference between the estimates based on bulk aearodynamical calculation and average volume and those with IMB non-steady state method was 1–2%, but for the smallest lake (Lehee) the difference was 5–7%.

Table 4. Evaporated fraction of water volume (f) in 5 lakes from 16 June to 18 August 2003 (Supplemental Fig. S5) estimated with different methods: (1) evaporation loss (E) estimated by bulk aerodynamic calculations (equation 4) and divided by lake average volumes (E/mean volume); (2) Lake Iso-Roine water level change-based estimate is also given (Water-level based) (data from Hertta database, https://www.syke.fi/avoindata, no data (nd) available for the other lakes); (3) isotope mass balance method (IMB no LEL adj) (IMB, see equation 6 in methods), and (4) with local evaporation line (LEL) adjustment (IMB LEL adj) (see Skrzypek et al. 2015).

Lake	f (E/mean volume)	f (Water-level based)	f (IMB no LEL adj)	f (IMB LEL adj)
Ormajärvi	0.02	nd	0.03	0.02
Suolijärvi	0.04	nd	0.05	0.03
Lehee	0.16	nd	0.10	0.07
Pyhäjärvi	0.02	nd	0.04	0.03
Iso-Roine	0.03	0.03	0.05	0.03

Discussion

Our study was conducted in a boreal landscape characterized by a high number of small- and medium-sized lakes (cf. Verpoorter et al. 2014) connected by streams and rivers to larger lakes downstream. Our results showed a rather high between-lake and within-lake variability of δ²H and δ¹⁸O values in a small geographical area. In the small headwater lakes, δ²H and δ¹⁸O values are closer to those in rain and groundwater, whereas in downstream lakes the enrichment of heavier isotopes occurs because of evaporation both at the lake surface and in the catchment (cf. Gibson and Edwards 2002, Gibson and Reid 2014). This finding is clearly seen in the epilimnetic δ¹⁸O values of the 10 study lakes measured in October after autumnal mixing, which were strongly related to lake surface area and proportion of inland waters of the catchment and showed ∼1‰ interannual difference between the dry and wet year (Fig. 6).

In addition to precipitation and runoff, hydraulic residence time of water both at the drainage basin and lake scale plays an important role for seasonal and interannual differences in isotope signatures of the lakes (Haig et al. 2021). Jonsson et al. (2009) showed that in small (area <1 km²) sub-Arctic lakes, seasonal variation in the isotopic composition of water can be clearly observed in lakes with short residence times (<6 months), whereas lakes with a longer residence time retain an isotopic signal closer to that of annual mean precipitation. Similarly, in our study lakes the greatest seasonal variation in δ²H and δ¹⁸O values occurred in the 2 smaller lakes (Lehee and Suolijärvi, area 1 and 2 km², respectively), with residence times <6 months. However, the impact of evaporation for seasonal isotopic variation was evident also in our largest lakes (area 6.5–31.5 km²) with longer residence times in summer 2003, when evaporation exceeded precipitation. In boreal lakes, the role of evaporation is more distinct than that observed in sub-Arctic lakes, in which the seasonal variation in δ²H and δ¹⁸O of lake waters is mainly influenced by winter precipitation (spring thaw) and summer precipitation (Jonsson et al. 2009). The evaporation estimates with bulk aerodynamic method in our study lakes (405–544 mm during May–Oct) were close to those obtained by Hyvärinen et al. (1995) and Stenberg (2007) for Lake Pääjärvi, a deep mesotrophic lake close to our study area (see review by Strachan et al. 2016 for comparison).

Although we did not have direct hydrological measurements in the lake chain, we could indirectly estimate water balances using discharge, precipitation, relative humidity, and temperature measurements in the nearby region and using IMB calculations. Ideally, IMB estimates are based on data of δ²H and δ¹⁸O values of all water sources and their progressive enrichment in steady state conditions due to evaporation downstream the landscape (Gibson et al. 2010, Gibson and Reid 2014). However, data for all water sources are not always available, and steady-state conditions do not always prevail, causing uncertainties in the results, especially during the warm and rather dry ice-free period (May–Oct) 2003, showing variable E/I estimates for the lakes. The IMB steady-state model produced unrealistic negative estimates for 2 lakes, although both evaporation at the lake surfaces and at least minor inflow occurred. During the warm and dry period when the lakes were strongly stratified from mid-June to mid-August 2003, however, the IMB nonsteady state model yielded similar results for the fraction evaporated volume of water from 5 lakes, especially the large lakes (Ormajärvi, Pyhäjärvi, Iso-Roine), compared with those obtained by using bulk aerodynamic evaporation estimation and water-level decrease (Iso-Roine) versus lake volume. Because the exact initial water volume in the lakes in June is not known and assumed as an average volume, the obtained results with different methods surprisingly support each other.

During the 2004 ice-free season, when evaporation was lower and precipitation higher than during the previous year, the E/I estimates based on bulk aerodynamic E and CA based inflow calculations and the ratio obtained by IMB steady-state method were closer each other; these independent results were strongly correlated. Interestingly, the IMB results using mean 1 m values for May–October and those sampled after water column mixing in October were similar and strongly correlated. Thus, in future water balance studies, sampling only in autumn (c.f. Bennett et al. 2008) may save effort and resources, although more detailed information, such as linear evaporative loss during June–August 2003, is then missed.

The IMB approach has proved useful for water balance estimations, especially in remote areas where direct hydrological measurements are infrequent or sparse, or even unavailable (e.g., Gibson and Edwards 2002, Bennett et al. 2008, Shi et al. 2017, Yapiyev et al. 2020). The greatest uncertainties in IMB estimates are in the highly sensitive calculations for isotopic values of ambient atmospheric moisture (Skrzypek et al. 2015). For this, we used average values (May–Oct) of δ²H and δ¹⁸O of precipitation in the study area. The observed strong seasonal variation in these values is a typical feature in northern Europe (Kortelainen and Karhu 2004, Gourcy et al. 2005, Jonsson et al. 2009), and therefore the E/I estimates based on stepwise interpolation of data collected every 2 weeks produced unrealistic estimates (data not shown). Bennett et al. (2008) proposed an adjustment factor to better estimate the isotope signature of ambient air moisture based on local precipitation and LEL, which Skrzypek et al. (2015) included in the Hydrocalculator software. In our study lakes, the E/I ratios calculated without this adjustment were closer to those obtained by bulk aerodynamic E- and CA-based E/I ratio, but we are unable to evaluate which values are closer to reality. The LEL slope in the study lakes (5.45) is typical of the lakes in a seasonal climate (Gibson et al. 2008).

Groundwater δ²H and δ¹⁸O values at our study area varied between −82.3–99.9‰ and 11.6–13.4‰, respectively, in line with those observed in southern Finland (Kortelainen and Karhu 2004). During winter, the effect of groundwater inflow is visible in many lakes, as indicated by holes in the ice near the shoreline (LA, pers. observ.). In Lake Pyhäjärvi, the lower δ¹⁸O values at 1 m depth than close to the bottom in winter may indicate groundwater inflow. Midwinter under ice δ¹⁸O values in Lake Ormajärvi also were lower close to the bottom than at the surface, suggesting warmer groundwater flow into the lake. In addition, higher evaporation to surface water inflow ratio (E/I) in Lake Ormajärvi than in the 2 large lakes (Pyhäjärvi and Iso-Roine), which had similar daily enrichment (0.009‰) of the δ¹⁸O values during the dry period in 2003, suggests that, especially in Lake Ormajärvi, our water budget estimates did not include all water inputs, although long residence time presumably influenced the results. Lake Kyynärö (outliers in Fig. 6a) with low δ¹⁸O values locates at the subcatchment of Lake Ormajärvi, and in both lakes groundwater presumably originates from the nearby esker. In addition, the higher isotope values in Lake Ormajärvi than in the next downstream lake (Suolijärvi) suggest additional groundwater inflow into this lake. In a boreal esker landscape (Rokua, northern Finland, area 90 km²) with more intensive measurements of δ²H and δ¹⁸O in groundwater, streams, and lakes, Isokangas et al. (2015) estimated the proportion of groundwater in the total water inflow to the lakes to vary between 39% and 98%. However, the application of the same method in our boreal lake chain in a much larger and complex drainage area did not yield realistic estimates, presumably because our measurements were from the deepest points of the lakes, and many groundwater release areas are unknown. Water in seepage lakes fed by groundwater is typically clear and oligotrophic, but groundwater flow from peatland areas can also contribute to nutrient and organic carbon loading to lakes and rivers (Ala-aho et al. 2013, Palviainen et al. 2021). In our study area, the importance of groundwater flow in the whole water mass transport is presumably much lower than in Rokua esker area (c.f. Arvola et al. 2015).

In conclusion, the knowledge of water origin and mass transport are crucial for management operations in lakes and their catchments. Because of intensifying land use and climate change with warmer winters and shortening ice-cover period (cf. Veijalainen et al. 2012, Ahopelto et al. 2019), water protection by reducing external loading from the watershed (e.g., from agriculture and forestry) is crucial. Hydrological conditions and pathways play significant roles for loading, transport, and retention of key elements (carbon, nitrogen, and phosphorus) and their metabolic processing in different subsystems within a drainage basin. Measurements of water stable isotope values provide tools for these estimates, but the results are not always unambiguous and sensitive for air moisture data. Our study suggests that, in addition to surface water runoff measurements, more geohydrological studies should be conducted to estimate groundwater flow to lakes. The data for this study were collected during years when the lakes were ice-covered and when summer drought and heavy precipitation periods differed from long-term average conditions. Thus, these results can provide valuable background reference information for future studies in warming climate and a shortening ice-cover period in the boreal landscape.

Acknowledgements

We thank Riitta Ilola and Jaakko Vainionpää, the laboratory technicians, and other staff members who helped us in the field and laboratory of Lammi Biologial Station, as well as the Chronology laboratory, University of Helsinki, for the facilities and support.

Disclosure statement

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Additional information

Funding

The study was part of the NUTRIBA project supported by the Academy of Finland within the research programme BIREME [Grant number202419].