Isotope Composition of Precipitation, Groundwater, and Surface and Lake Waters from the Plitvice Lakes, Croatia

: The application of tritium, 2 H, and 18 O in the characterization of the precipitation, groundwater, and surface and lake water of the Plitvice Lakes (PL), Croatia, over the 1979–2019 period is presented. An increase in the mean annual air temperature of 0.06 ◦ C / year and in the annual precipitation amount of 10 mm / year is observed. The good correlation of the tritium activity concentration in the PL and Zagreb precipitation implies that the tritium data for Zagreb are applicable for the study of the PL area. The best local meteoric water line at PL was obtained by the reduced major axis regression (RMA) and precipitation-weighted ordinary least squares regression (PWLSR) approaches: δ 2 H PWLSR = (7.97 ± 0.12) δ 18 O + (13.8 ± 1.3). The higher deuterium excess at PL (14.0 ± 2.2 % (cid:24) ) than that at Zagreb reﬂects the higher altitude and inﬂuence of the Mediterranean precipitation. The δ 2 H in precipitation ranges from − 132.4% (cid:24) to − 22.3% (cid:24) and δ 18 O from − 18.3 % (cid:24) to − 4.1% (cid:24) . The much narrower ranges in the groundwater ( < 1% (cid:24) in δ 18 O, < 10% (cid:24) in δ 2 H) indicate the good mixing of waters in aquifers and short mean residence times. The higher average δ 2 H in all three karst springs observed after 2003 can be attributed to the increase in the mean air temperature. The mean δ 2 H and δ 18 O values in the surface and lake water increase downstream due to the evaporation of surface waters. There is no signiﬁcant di ﬀ erence between the surface water line and the lake water line (2011–2014). The stable isotope composition of the surface and lake waters reacts to extreme hydrological conditions.


Introduction
Karst is a special type of landscape that is formed by the process of karstification-i.e., the dissolution of soluble carbonate rocks, mostly limestone and dolomite. Karst channels, conduits, and fissures store relatively large quantities of groundwater, and such karst aquifers are capable of providing large supplies of water for human consumption. The precipitation water quickly infiltrates underground, creating a system of interconnected flow paths, and eventually re-appears at the surface as springs. Karst aquifers, especially in areas with a high permeability, can be very vulnerable to contamination and can enable the fast transport of contaminants through the aquifers, which can result in the degradation of water quality [1][2][3][4]. The assessment of the impact of human activities and recent climate changes on karst waters has to be properly considered [5]. A good understanding The Plitvice Lakes (PL) area presents a unique system of 16 cascade flow-through lakes that are fed by three main springs (Crna Rijeka River, Bijela Rijeka River, Plitvica River) and outflow to the Korana River. It is famous worldwide for its beauty and diversity. The area has been protected as a part of the national park since 1949, and since 1979 it has been included in the United Nations Educational, Scientific, and Cultural Organisation (UNESCO) World Heritage List. The area is one of the most studied karst areas in Croatia, where various scientific studies have been performed since the beginning of 20th century. During the 1970s, a group from the Ruđer Bošković Institute joined

Site Description
The Plitvice Lakes area is located in a continental, mountainous part of west Croatia (Figure 1). The Plitvice Lakes are a system of 16 lakes developed on carbonate rocks [41] separated by tufa barriers and interconnected by waterfalls ( Figure 2). Geologically, they belong to the External Dinarides or Dinaric-coastal area, with characteristic carbonate sediments and karstic features. The altitude difference of the system from the springs to the Korana River is about 300 m (Figure 2b). However, the recharge area extends to higher mountains (the highest altitude of the national park is 1279 m a.s.l.), and the average recharge altitude is about 900 m a.s.l. [38]. Zagreb is also situated in a continental part of Croatia (to the north) ( Figure 1). It belongs to the Pannonian area characterized by its milder relief; predominantly magmatic, clastic, and metamorphic rocks; and well-developed stream grid [42]. Zadar is located at the Adriatic Sea coast (southern Croatia) (Figure 1), which belongs to the Dinaric coastal area.
The three selected locations of precipitation sampling differ in the altitude of the meteorological stations; Plitvice Lakes station is situated at 550 m a.s.l., Zagreb-Grič station at 165 m a.s.l., and Zadar at 5 m a.s.l. [19,23]. Climatologically, Zagreb and Plitvice Lakes belong to the Cfb climate class, characterized by a temperate climate without a dry season and with a warm summer, while Zadar belongs to the Cfa climate class, characterized by a temperate climate without a dry season but with a hot summer [43,44]. The annual precipitation amounts in Zagreb (measured at the Zagreb-Grič meteorological station) and Zadar are almost identical: 883 mm and 915 mm for the  period, and 884 mm and 882 mm for the 1991-2004 period, respectively [19]. The annual precipitation amount is significantly higher at the Plitvice Lakes, where it ranges between 1148 and 2113 mm in the 1986-2019 period [45]. The monthly precipitation at the Plitvice Lakes is distributed relatively uniformly throughout the year, with slight maxima observed in spring and autumn and minimum values in summer months. The annual temperatures range from 8.0 to 10.8 • C (1986-2019), with an average value of 9.2 ± 0.5 • C. January is the coldest month (0 ± 2.3 • C on average) and July is the warmest month (18.4 ± 1.1 • C) [45]. Snow falls between November and March; however, reduced snowfall is observed in recent times compared to the older data [33,45].
The water temperature at the springs is very stable throughout the year. Nevertheless, an increase in the spring temperatures is observed between the 1981-1986 and 2010-2014 periods in both Crna Rijeka spring (from 7.80 ± 0.15 • C to 8.04 ± 0.14 • C) and the Bijela Rijeka spring (from 7.46 ± 0.17 • C to 8.14 ± 0.53 • C) [33]. The temperature of the surface water downstream ranges from near zero in winter to 23 • C in summer, being on the average about 12 • C [33,46]. The lakes can be frozen in the winter months.  Table 1.   Table 1.

Sampling
The monthly precipitation was collected between 1978 and 1984 [36] (not continuously) and from 2003 to 2006 [11,37] at the meteorological station at the Plitvice Lakes (altitude 550 m, International Atomic Energy Agency (IAEA) and World Meteorological Organization (WMO) station code 1432501) [23].
Groundwater samples were collected at the three main karst springs of the system: Crna Rijeka River and Bijela Rijeka River (in the south), and the Plitvica River that joins the lake water after the Great waterfall before the location of Korana River-Sastavci ( Figure 2, Table 1).
Surface water was collected as grab samples at 8 locations along the water course for the length of~10 km from Matica (the main stream feeding the lakes) to the Korana River (the outflow from the lake system) ( Figure 2, Table 1). Lake water was collected at 4 sediment traps in different lakes: IRB1 and IRB2 in Lake Kozjak at water depths of 6 m and 8-10 m, respectively; IRB3 in Lake Prošćansko at 6 m; and IRB4 in Lake Gradinsko at 2 m water depths (note: the water from sediment traps will be called "lake water" in the further text). At the same time, the surface water was collected at nearby locations.
All the sampling locations are listed in Table 1, with the location code and the type of water body sampled. The amount of tritium activity concentration data, stable isotope data, and period of sampling is given in Table S1.

Meteorological and Hydrological Data
The meteorological data consisted of the monthly precipitation amount (P) and the average monthly air temperature (T). The data were obtained on request from the Croatian Meteorological and Hydrological Service (CMHS) [45]. The meteorological records for the Plitvice Lakes, Croatia, exist for three distinctive periods: 1986-1990, 1996-2011, and 2015-2019. The data are not complete for all years, and in further analyses only the years with all 12 months of data will be used. The minimal and maximal monthly values within a year were identified, and the mean annual temperature and total annual precipitation amount were determined. CMHS [45] also provided data on the flow rates measured at different points in the Plitvice Lakes area for the period since 1982, except for the 1991-2001 period.

Measurement
Details on the measurement techniques and their changes were described in [23], and here we give only a short overview. The stable isotope composition of water samples for the period up to 2003 was measured on a Varian MAT 250 dual inlet isotope ratio mass spectrometer (IRMS) at the Jožef Stefan Institute in Ljubljana [23,47]. The measurement precision of duplicates was better than ± 0.1% for δ 18 O and ± 1% for δ 2 H. The δ 2 H and δ 18 O in the surface and lake water samples in the period 2012-2014 were measured at the JOANNEUM, Graz, Austria. The oxygen isotopic composition was determined on a dual-inlet Finnigan DELTAplus by means of the fully automated equilibration technique, and the isotopic composition of hydrogen was determined on a continuous flow Finnigan DELTAplus XP mass spectrometer with a HEKAtech high-temperature oven by the reduction of water over hot chromium [48]. All the measurements were carried out together with laboratory standards that were calibrated periodically against international standards, as recommended by the IAEA. The measurement precision was better than ± 0.1% for δ 18 O and ± 1% for δ 2 H [47]. The stable isotope composition of recent samples was determined at the Laboratory for Spectroscopy of the Faculty of Mining, Geology, and Petroleum Engineering, University of Zagreb, with a Liquid Water Isotope Analyzer (LWIA-45-EP, Los Gatos Research), and the official LGR working standards were used. The data were analyzed by the Laboratory Information Management System (LIMS) [49].
The tritium activity concentration (A) in all the samples was determined at the Ruder Bošković Institute in Zagreb, except for some data for the 2003-2006 period [10]. The results are expressed in tritium units (1 TU = 0.118 Bq l −1 ) [20,22], which represent one 3 H atom in 10 18 atoms of hydrogen. The gas proportional counting technique (GPC) was used up to 2009 [34,50]. The detection limit was 2.5 TU, and the measurement uncertainty was between 2 and 5 TU, depending on the activity concentration. Since 2008, the technique of the electrolytical enrichment of samples (LSC-EE) and measurement by an ultra-low-level liquid scintillation counter Quantulus 1220 was used [51]. The detection limit obtained by the LSC-EE technique was 0.5 TU, and the measurement uncertainty was between 0.5 and 3 TU [51].

Data Analysis
The results of the stable isotope analyses are reported as δ-values-i.e., the relative difference in isotope ratios of the sample and the standard [52]: where R Sample and R Reference stand for the isotope ratio (R = 2 H/ 1 H and R = 18 O/ 16 O) in the sample and the reference material (standard), respectively. The δ-values are dimensionless and small, and therefore they are expressed in per mill (% ) [21,[52][53][54]. The international standard VSMOW (the Vienna Standard Mean Ocean Water) is used [23,54,55]. The δ 2 H and δ 18 O isotopic compositions of meteoric waters (precipitation and atmospheric water vapor) are strongly correlated, and the relation in Equation (2) is referred to as the global meteoric water line (GMWL) [53,[55][56][57].
The GMWL describes the general relation between δ 2 H and δ 18 O on a global scale reasonably well. However, for applications in hydrogeological studies, regional local meteoric water lines (LMWLs), either long-term or for certain shorter periods, can be more appropriate [23,58,59]. Generally, a LMWL has the form δ 2 H = a δ 18 O + b, where a is the slope and b is the intercept. LMWLs can differ from the GMWL in terms of both the slope and intercept values, depending on the conditions for forming a local water source [21,[58][59][60].
The deuterium excess (d-excess, or d) was calculated from the paired monthly data according to Equation (3) [56]: This can be related to the meteorological conditions in the source region from which the water vapor is obtained [20,21,[58][59][60]. Autumn and winter precipitation originating from the Mediterranean Sea is characterized by distinctly higher d-excess values (d > 18% ) than precipitation coming from the Atlantic (d~10% ), reflecting the specific source conditions during water vapor formation [19,59,60].
Correlations between various data points were obtained as ordinary least squares regressions using standard commercial software. Pearson's coefficient r is given, as are the number of data pairs n and the p-values describing the statistical significance of the correlations. The data taken from the literature usually have the adjacent R 2 value reported.
In the special case of calculating the aand b-values of LMWLs, different methods were applied: ordinary least squares regression (OLSR), reduced major axis regression (RMA), and major axis least squares regression (MA, or the orthogonal regression) [61,62]. We calculated precipitation-weighted regressions (PWLSR, PWRMA, and PWMA) [61][62][63], which took into account the precipitation amount in a particular month. The local meteoric water lines are defined as LMWL OLSR , LMWL RMA , LMWL MA , LMWL PWLSR , LMWL PWRMA , and LMWL PWMA . For calculating the regressions from data sets of at least 36 continuous monthly records, we used the Local Meteoric Water Line Freeware [64]. The software also calculated an average of the root mean square sum of squared errors (rmSSEav), which is a relative error that allows for a comparison of different methods; the closer the value of rmSSEav is to 1.0, the better the regression method is for that set of data [61].

Results and Discussion
In this section, we present an analysis of the climatological parameters (temperature T and precipitation amount P) in the area of the Plitvice Lakes, followed by the isotope composition (tritium activity concentration, δ 2 H and δ 18 O values) of precipitation, groundwater, and finally the surface and lake waters.

Temperature and Precipitation Amount
The monthly mean temperatures in three periods when data are available-1986-1990, 1996-2011, and 2015-2019-show an increase in temperature in all months (Figure 3a). The increase is statistically significant at 5% (p < 0.05) in June, July, August, and November, and at 10% (p < 0.10) in April. When the mean annual temperatures in the 1986-2019 period are compared (Figure 3b) (only years with data for all 12 months), an increase with the slope of 0.06 ± 0.01 • C per year is observed, and n = 16, r = 0.85, p < 0.05. The mean annual temperature at the Zagreb-Grič station in the 1976-2018 period showed a significant increase, with a slope of 0.071 ± 0.008 • C per year (r = 0.82, p < 0.05) and a faster increase in the maximal annual temperatures (0.09 ± 0.02 • C per year, r = 0.69) [23], which is in accordance with the observation at the Plitvice Lakes, although the set of meteorological data at the Plitvice Lakes is not complete. The mean temperatures in the three periods with available data ( Table 2) show a constant increase. (Table 2) show a constant increase.  The annual increases in temperature at both the Zagreb and PL stations can be compared with the global temperature increase and the temperature increase in the city of Ljubljana in Slovenia. Ljubljana shows a distinctive air warming trend, particularly in the period 1979-2008, of 0.06 °C per year [65]. Although this period does not include the hottest years on record globally (2014-2018), with 2016 being the hottest year [66], the value of the increase is consistent with the values obtained for Zagreb [23] and the Plitvice Lakes. Moreover, all the data are higher than the globally observed temperature change of 0.018 °C per year in the 1980-2020 period [67]. A recent study [68] indicated that the cities in East Europe, including Ljubljana, will experience air warming at a higher extent than is observed globally.  The monthly mean precipitation amounts show no significant seasonal variations (Figure 4a), although a slight minimum is observed in the summer months, and slightly higher values in spring and autumn. Higher fluctuations within a month are observed in the most recent period (2015-2019) compared to the older ones ( Figure 4a). The annual precipitation amounts increase with a slope of 10 ± 7 mm per year, r = 0.31 ( Figure 4b, Table 2). A wider range in the monthly and annual precipitation amount values has been observed ( Figure 4b, Table 2), as observed also in the data for Zagreb [23]. The long-term (1976-2018) annual precipitation amount at Zagreb showed a slight, statistically not significant, increase of 1.4 ± 1.7 mm per year. However, the most prominent characteristics of the data were the higher deviations of the annual values in the period after 2000 from the mean value for the whole period [23].  The annual increases in temperature at both the Zagreb and PL stations can be compared with the global temperature increase and the temperature increase in the city of Ljubljana in Slovenia. Ljubljana shows a distinctive air warming trend, particularly in the period 1979-2008, of 0.06 • C per year [65]. Although this period does not include the hottest years on record globally (2014-2018), with 2016 being the hottest year [66], the value of the increase is consistent with the values obtained for Zagreb [23] and the Plitvice Lakes. Moreover, all the data are higher than the globally observed temperature change of 0.018 • C per year in the 1980-2020 period [67]. A recent study [68] indicated that the cities in East Europe, including Ljubljana, will experience air warming at a higher extent than is observed globally.
The monthly mean precipitation amounts show no significant seasonal variations (Figure 4a), although a slight minimum is observed in the summer months, and slightly higher values in spring and autumn. Higher fluctuations within a month are observed in the most recent period (2015-2019) compared to the older ones ( Figure 4a). The annual precipitation amounts increase with a slope of 10 ± 7 mm per year, r = 0.31 ( Figure 4b, Table 2). A wider range in the monthly and annual precipitation amount values has been observed ( Figure 4b, Table 2), as observed also in the data for Zagreb [23]. The long-term (1976-2018) annual precipitation amount at Zagreb showed a slight, statistically not significant, increase of 1.4 ± 1.7 mm per year. However, the most prominent characteristics of the data were the higher deviations of the annual values in the period after 2000 from the mean value for the whole period [23].

Isotope Composition of Precipitation
The record of the tritium activity concentration (A) in the precipitation for Plitvice Lakes (Table  S2) is not as complete as the one for Zagreb precipitation [23]  It is obvious ( Figure 5) that the tritium activity concentration in the precipitation at Plitvice Lakes follows the trend of the precipitation in Zagreb. The linear correlation between the two sets revealed a line with a slope of (1.02 ± 0.04), n = 62, r = 0.96, and the paired t-test showed that at the 0.05 level, the two sets of data are not significantly different.

Isotope Composition of Precipitation
The record of the tritium activity concentration (A) in the precipitation for Plitvice Lakes (Table S2) is not as complete as the one for Zagreb precipitation [23]  It is obvious ( Figure 5) that the tritium activity concentration in the precipitation at Plitvice Lakes follows the trend of the precipitation in Zagreb. The linear correlation between the two sets revealed a line with a slope of (1.02 ± 0.04), n = 62, r = 0.96, and the paired t-test showed that at the 0.05 level, the two sets of data are not significantly different.
The  Table 2) can account for about 1% difference in the δ 18 O values, taking into account the temperature gradient of 0.331% per • C [23].
A difference in altitude of 385 m and the altitude effect of 0.28% per 100 m altitude difference [19] would account for a further 1.1% , giving a total difference of 2.1% in δ 18 O, which is in good agreement with the observed difference, taking into account limited number of isotope data for the Plitvice Lakes.  Table 2) can account for about 1‰ difference in the δ 18 O values, taking into account the temperature gradient of 0.331‰ per °C [23]. A difference in altitude of 385 m and the altitude effect of 0.28‰ per 100 m altitude difference [19] would account for a further 1.1‰, giving a total difference of 2.1‰ in δ 18 O, which is in good agreement with the observed difference, taking into account limited number of isotope data for the Plitvice Lakes.  Table 3. Values of slopes (a ± σ a ) and intercepts (b ± σ b ) for LMWLs obtained by different regression methods for precipitation at the Plitvice Lakes, Zadar, and Zagreb. OLSR: ordinary least squares regression; RMA: reduced major axis regression; MA: major axis least squares regression; PWLSR, PWRMA, and PWMA: precipitation-weighted respective regressions; n-number of data points included; rmSSEav-average of the root mean square sum of squared errors [60]. In bold: approaches that describe LMWL the best.

Location and Period
Method   In order to determine the relation between the δ 2 H and δ 18 O values in the precipitation of Plitvice Lakes-i.e., the LMWL of the form δ 2 H = (a ± σa) δ 18 O + (b ± σb)-we performed linear regression by different approaches (Table 3). For comparison, the analysis for the Zagreb precipitation is taken from [23], while the data for Zadar [19] were analyzed here also by applying the same approaches ( Table 3). The analysis of the stable isotope composition of the Plitvice Lakes precipitation is based on 36 data pairs from 2003-2006 (Table 1, Table S2). The LMWLs obtained by various correlation methods for all three precipitation stations (PL, Zg, and Zd) are shown in Table 3. The best results for the Plitvice Lakes and Zadar precipitation were obtained by the RMA and PWLSR approaches (marked in bold in Table 3), while the Zagreb precipitation data were best described by the RMA and PWMA approaches (marked in bold in Table 3). The LMWLs obtained by the best approaches are shown in Figure 6. The LMWLRMA and LMWLPWLSR for PL have slopes close to slope 8 of the Global Meteoric Water Line (GMWL), but slightly higher intercepts. The slope of the Zd LMWLs is lower, as well as the intercepts, caused by evaporation during the summer months [19]. In general, the Plitvice Lakes LMWL lines lie above the GMWL and LMWLs for Zagreb, and are close to the LMWLs for Zadar ( Figure 6). Several precipitation sampling campaigns of precipitation related to the cave drip water studies in the karst area of Croatia revealed similar results for LMWLs: sites along the coast at low altitudes gave values of slope (6.6-6.8) and intercept (3.8-6.7) similar to those obtained for Zadar, while continental stations gave slopes (7.1 to 7.8) and intercepts (2.3 to 14.5) closer to the ones of Zagreb and the Plitvice Lakes [25]. Table 3. Values of slopes (a ± σa) and intercepts (b ± σb) for LMWLs obtained by different regression methods for precipitation at the Plitvice Lakes, Zadar, and Zagreb. OLSR: ordinary least squares regression; RMA: reduced major axis regression; MA: major axis least squares regression; PWLSR, PWRMA, and PWMA: precipitation-weighted respective regressions; n-number of data points included; rmSSEav-average of the root mean square sum of squared errors [60]. In bold: approaches that describe LMWL the best.  (Table 3). For comparison, the analysis for the Zagreb precipitation is taken from [23], while the data for Zadar [19] were analyzed here also by applying the same approaches ( Table 3). The analysis of the stable isotope composition of the Plitvice Lakes precipitation is based on 36 data pairs from 2003-2006 (Table 1, Table S2). The LMWLs obtained by various correlation methods for all three precipitation stations (PL, Zg, and Zd) are shown in Table 3. The best results for the Plitvice Lakes and Zadar precipitation were obtained by the RMA and PWLSR approaches (marked in bold in Table 3), while the Zagreb precipitation data were best described by the RMA and PWMA approaches (marked in bold in Table 3). The LMWLs obtained by the best approaches are shown in Figure 6. The LMWL RMA and LMWL PWLSR for PL have slopes close to slope 8 of the Global Meteoric Water Line (GMWL), but slightly higher intercepts. The slope of the Zd LMWLs is lower, as well as the intercepts, caused by evaporation during the summer months [19]. In general, the Plitvice Lakes LMWL lines lie above the GMWL and LMWLs for Zagreb, and are close to the LMWLs for Zadar ( Figure 6). Several precipitation sampling campaigns of precipitation related to the cave drip water studies in the karst area of Croatia revealed similar results for LMWLs: sites along the coast at low altitudes gave values of slope (6.6-6.8) and intercept (3.8-6.7) similar to those obtained for Zadar, while continental stations gave slopes (7.1 to 7.8) and intercepts (2.3 to 14.5) closer to the ones of Zagreb and the Plitvice Lakes [25].
The deuterium-excess monthly values at the Plitvice Lakes range from 7.7% (January 2004) to 17.9% (January 2005), with the mean value of 14.0 ± 2.2% . These values are higher than those for the Zagreb precipitation, where the mean annual d-excess values ranged between 2.3% and 10.7% , and in the 2003-2006 period the mean value was 8.8 ± 0.8% [23]. The higher d value may be explained by the observed increase in the d-excess value with the altitude [2,19], but also by the more intense influence of the Mediterranean air masses in the area of Plitvice lakes [13,19], which is shown here by the relative position of the LMWL compared to the GMWL ( Figure 6). Similar high d-excess values (between 11.0 and 18.6) were found for several precipitation sampling locations in the continental karst area at altitudes between 300 and 1550 m [25].
The deuterium excess values for the Zagreb precipitation were higher in autumn than in spring. [19,23,69]. The higher d-excess in autumn indicates a higher influence of the Mediterranean air masses in these months. A similar pattern of the monthly d-excess distribution was observed for the Plitvice Lakes station (Figure 7), with higher d-excess monthly values observed in the autumn-winter precipitation (September-December, 15.2 ± 0.4% ). The lowest mean monthly d-excess values were observed in the summer (July-August, 12.5 ± 1.0% ).

Groundwater
The tritium activity concentration in the three main springs of the Plitvice Lakes ( Figure 8) shows a general decrease over the sampling period, but the seasonal variations are much smaller than those for precipitation ( Figure 5), and its values fluctuate around the mean annual A values for precipitation. The A values in the Bijela Rijeka springs were always higher than in the other two springs. It should be mentioned here that the δ 18 O values from the 1983-1984 period ranged from −7.2% to −13.1% , the δ 2 H from −46.0% to −94.4% , and the d-excess from 8.9% to 11.5% , based on only five data points for the monthly isotope composition of the precipitation at the Plitvice Lakes [7,8].

Groundwater
The tritium activity concentration in the three main springs of the Plitvice Lakes ( Figure 8) shows a general decrease over the sampling period, but the seasonal variations are much smaller than those for precipitation ( . In addition, the relative standard deviations of 7.8% in 1984 and 1.3% in 2015 in A for BR were smaller in both periods than those for CR (17% and 9%, respectively).
The tritium activity concentration in the three main springs of the Plitvice Lakes ( Figure 8) shows a general decrease over the sampling period, but the seasonal variations are much smaller than those for precipitation ( Figure 5), and its values fluctuate around the mean annual A values for precipitation. The A values in the Bijela Rijeka springs were always higher than in the other two springs. The range in the δ 18 O in precipitation was 14.2‰ and in δ 2 H it was 110.1‰ (Table S2). However, the ranges in the δ values in groundwater (at springs) were much lower, at less than 1‰ in δ 18 O and less than 10‰ in δ 2 H ( Figure 10  The range in the δ 18 O in precipitation was 14.2% and in δ 2 H it was 110.1% (Table S2). However, the ranges in the δ values in groundwater (at springs) were much lower, at less than 1% in δ 18 O and less than 10% in δ 2 H ( Figure 10  Although the isotope data presented above and the physico-chemical and carbon isotope data from previous studies [7,33,46] indicate that the three springs originated from different hydrogeological units, and in principle different groundwater water lines should be obtained, the δ 2 H and δ 18 O data in the three springs cluster together along the LMWL line for the Plitvice Lakes, and the Groundwater Water Line (GWL) is determined for all the springs ( Figure 11, Table 4): The GWL (Equation (4)) has both a lower slope and a lower intercept than the LMWL (Table 3), which may be a consequence of water evaporation during its flow in karst aquifers, but also due to a very narrow range of the δ 2 H values for BR. The regional groundwater line (RGWL) for the Croatian karst area, based on 340 karst springs, revealed a slope of 7.0 and an intercept of 4.2 (Table 4) [30]. Table 4. Values of the slopes (a ± σa) and intercepts (b ± σb) for the groundwater line (GWL), surface water lines (SWL), and lake water lines (LWL).  Although the isotope data presented above and the physico-chemical and carbon isotope data from previous studies [7,33,46] indicate that the three springs originated from different hydrogeological units, and in principle different groundwater water lines should be obtained, the δ 2 H and δ 18 O data in the three springs cluster together along the LMWL line for the Plitvice Lakes, and the Groundwater Water Line (GWL) is determined for all the springs ( Figure 11, Table 4):  (Table 3).

Type of Line
The physico-chemical parameters in the spring waters (temperature, pH, ion concentrations) were also very constant throughout year, and all these together indicated a good mixing of water, uniformly recharged during the year in ground water karst aquifers, and the mean residence time of Figure 11. Groundwater water line (GWL) for the three springs, all data after 2000, GWL is given by Equation (4). Comparison with the LMWL PWLSR for the Plitvice Lakes (Table 3). Table 4. Values of the slopes (a ± σ a ) and intercepts (b ± σ b ) for the groundwater line (GWL), surface water lines (SWL), and lake water lines (LWL). The GWL (Equation (4)) has both a lower slope and a lower intercept than the LMWL (Table 3), which may be a consequence of water evaporation during its flow in karst aquifers, but also due to a very narrow range of the δ 2 H values for BR. The regional groundwater line (RGWL) for the Croatian karst area, based on 340 karst springs, revealed a slope of 7.0 and an intercept of 4.2 (Table 4) [30].

Type of Line
The physico-chemical parameters in the spring waters (temperature, pH, ion concentrations) were also very constant throughout year, and all these together indicated a good mixing of water, uniformly recharged during the year in ground water karst aquifers, and the mean residence time of several years [7,33,46]. The conclusion was confirmed by applying the tritium activity concentrations in groundwater/springs for the determination of the Mean Residence Time (MRT) by the exponential model, which supposed a complete mixing of a new recharge with the already existing water in the aquifer [6][7][8]70]. The measured A values for the Zagreb precipitation were used as the input data. The MRT values of different karst springs were very short [8], ranging between 1 and 4 years on average. Among the three springs in the Plitvice Lakes area, the Crna Rijeka spring (CR) had the shortest MRT (2 years) and the Bijela Rijeka spring (BR) had the longest of about 4 years. All these data corroborate the previously obtained isotope data: the longer the MRT, the smaller the range of the isotope composition of groundwater ( Figure 9).
The extreme climatological characteristics in 1983 and 1984 that enabled the determination of the relative contributions of base-flow and precipitation in the three karst springs are of special interest here. Summer and fall/autumn 1983 were extremely warm and dry, so the aquifers were not completely recharged, and as a result the older water with higher tritium activity concentrations appeared at the springs (Figure 12a). In the following autumn, abundant precipitation was recorded, as well as a high amount of snow in the winter months (February and March 1984), and these caused lower tritium activities in springs in the spring of 1984 (Figure 12a). There are no available P and T data at the Plitvice Lakes in 1983 and 1984. However, the data for Zagreb could help describe the climatological conditions in these years. The mean annual temperature in Zagreb in 1983 was 12.1 • C, which is higher than that in 1982 (11.7 • C) and in 1984 (10.9 • C), and also higher than the average temperature for the 1980-1985 period of 11.2 ± 0.9 • C [23]. The precipitation amount in 1983 was 755 mm, which is lower than that in 1982 (805 mm) and in 1984 (897 mm), and also lower than the mean P in the 1980-1985 period (843 ± 67 mm) [23]. The flow rates at all three springs in autumn 1983 were very low (the lowest flow rates in November 1983: 0.36, 0.018, and 0.12 m 3 s −1 for CR, BR, and PR, respectively), much lower than the average flow rates in 1983 (1.327, 0.387, and 0.407 m 3 s −1 for CR, BR, and PR, respectively) or the long-term averages (1.41, 0.466, and 0.668 m 3 s −1 ). At all three springs, the highest flow rates were measured in April 1984 (9.12, 1.81, and 13.3 m 3 s −1 for CR, BR, and PR, respectively) [45] (Figure 12b).  The exponential model was applied for MRT calculation under the conditions of the highest and lowest A values in springs, although strictly the exponential model is not valid for non-uniform recharge. The calculated MRT from the higher A values in summer was doubled, for CR to 4 years and for BR to 8 years, in comparison with the previously determined 2 and 4 years, respectively. A significant decline in the tritium activity concentrations in springs was observed after abundant autumn and winter precipitation with lower tritium activity concentrations [7,8,35]. The lowest A values were observed at all three springs in April 1984 after intensive snow melting, and the corresponding MRTs were lower, at~1 year and~2 years for CR and BR, respectively. The lowest δ 2 H values in spring in April 1984 ( Figure 10) reflected the high contribution of winter precipitation. Knowing the A of precipitation, and the A in springs at the highest values (during the dry period) and at the lowest ones (April 1984), a fraction of precipitation input in the groundwater was estimated. The fraction of new (precipitation) water in April 1984 was the highest in CR, at about 90%, and was about 60% in PR, while BR consisted of approximately equal proportions of old groundwater and new input (50% each).
A similar conclusion about the mixing of quick and slow components in the springs of the Plitvice Lakes area was obtained by multiple tracers [10,11]. To obtain the MRTs, a multi-tracer lumped parameter modelling approach was applied using the time series of stable isotopes ( 2 H and 18 O) and tritium with a monthly resolution, as well as noble gases and CFCs. The average MRT of most springs was less than 5 years [10,11]. Two components of the groundwater flow were distinguished: the quick flow in the conduit network with an MRT of up to 0.5 years, and the slow component of the matrix of a fissured-porous aquifer with an MRT of up to 28 years [10,11].

Surface and Lake Waters
A relatively large set of δ 2 H data in surface waters in the period 2003-2005 can be found in [10], but the number of δ 18 O data is significantly lower (Table S1). The data show seasonal fluctuations, with maxima in late summer and minima in late winter (Figure 13a). However, the amplitude in δ 2 H of less than 10% is much smaller than the amplitude of the δ 2 H values in precipitation (about 110% , Figure 6), and reflects the isotope composition of groundwaters ( Figure 6). The maxima are not observed at all locations at the same time due to the retention of water in the two large lakes (Lake Prošćansko and Lake Kozjak). The two subsequent locations, LB before Lake Kozjak and KzB at the end of Lake Kozjak, show a very similar seasonal pattern, with a delay of approximately two months. The data are almost identical at locations KzB and KoB, although four smaller lakes are situated in between, indicating a fast flow of water through the small lakes. The average values increase in the downstream direction (Figure 13b) due to the evaporation of the surface waters, resulting in a Surface Water Line (SWL) with a lower slope (5.0 ± 0.6) and intercept (−17 ± 6) than the PL LMWL PWLSR line (Table 4).
Water 2020, 12, x FOR PEER REVIEW 18 of 27 and for BR to 8 years, in comparison with the previously determined 2 and 4 years, respectively. A significant decline in the tritium activity concentrations in springs was observed after abundant autumn and winter precipitation with lower tritium activity concentrations [7,8,35]. The lowest A values were observed at all three springs in April 1984 after intensive snow melting, and the corresponding MRTs were lower, at ~1 year and ~2 years for CR and BR, respectively. The lowest δ 2 H values in spring in April 1984 ( Figure 10) reflected the high contribution of winter precipitation.
Knowing the A of precipitation, and the A in springs at the highest values (during the dry period) and at the lowest ones (April 1984), a fraction of precipitation input in the groundwater was estimated. The fraction of new (precipitation) water in April 1984 was the highest in CR, at about 90%, and was about 60% in PR, while BR consisted of approximately equal proportions of old groundwater and new input (50% each). A similar conclusion about the mixing of quick and slow components in the springs of the Plitvice Lakes area was obtained by multiple tracers [10,11]. To obtain the MRTs, a multi-tracer lumped parameter modelling approach was applied using the time series of stable isotopes ( 2 H and 18 O) and tritium with a monthly resolution, as well as noble gases and CFCs. The average MRT of most springs was less than 5 years [10,11]. Two components of the groundwater flow were distinguished: the quick flow in the conduit network with an MRT of up to 0.5 years, and the slow component of the matrix of a fissured-porous aquifer with an MRT of up to 28 years [10,11].

Surface and Lake Waters
A relatively large set of δ 2 H data in surface waters in the period 2003-2005 can be found in [10], but the number of δ 18 O data is significantly lower (Table S1). The data show seasonal fluctuations, with maxima in late summer and minima in late winter (Figure 13a). However, the amplitude in δ 2 H of less than 10‰ is much smaller than the amplitude of the δ 2 H values in precipitation (about 110‰, Figure 6), and reflects the isotope composition of groundwaters ( Figure 6). The maxima are not observed at all locations at the same time due to the retention of water in the two large lakes (Lake Prošćansko and Lake Kozjak). The two subsequent locations, LB before Lake Kozjak and KzB at the end of Lake Kozjak, show a very similar seasonal pattern, with a delay of approximately two months. The data are almost identical at locations KzB and KoB, although four smaller lakes are situated in between, indicating a fast flow of water through the small lakes. The average values increase in the downstream direction (Figure 13b) due to the evaporation of the surface waters, resulting in a Surface Water Line (SWL) with a lower slope (5.0 ± 0.6) and intercept (−17 ± 6) than the PL LMWLPWLSR line (Table 4).   The more recent sampling campaign (between 2011 and 2014) included also lake water from sediment traps at certain water depths (Table 1, Figure 14a, Figure S1). No difference between lake water at selected depths and surface waters at the closest sampling sites is observed when the mean values are concerned, as can be seen from basic statistics data presented in Table S3. A slight increase in both the mean δ 2 H and δ 18 O values and their seasonal variations is observed for locations along the water course. The δ 18 O increased from −10.7 ± 0.1% at Matica to −10.3 ± 0.2% at the Korana River-Sastavci (KoS) ( Figure S2). Similarly, the δ 2 H increased from −71.4 ± 1.1% at Matica to −69.0 ± 1.7% at KoS (Figure 14b). Both the δ 2 H and δ 18 O are slightly lower at the final location in this study, KoB, probably due to some local groundwater input along the Korana River course. There is no significant difference between the isotope composition of the lake water and the corresponding surface water. The more recent sampling campaign (between 2011 and 2014) included also lake water from sediment traps at certain water depths (Table 1, Figure 14a, Figure S1). No difference between lake water at selected depths and surface waters at the closest sampling sites is observed when the mean values are concerned, as can be seen from basic statistics data presented in Table S3. A slight increase in both the mean δ 2 H and δ 18 O values and their seasonal variations is observed for locations along the water course. The δ 18 O increased from −10.7 ± 0.1‰ at Matica to −10.3 ± 0.2‰ at the Korana River-Sastavci (KoS) ( Figure S2). Similarly, the δ 2 H increased from −71.4 ± 1.1‰ at Matica to −69.0 ± 1.7‰ at KoS (Figure 14b). Both the δ 2 H and δ 18 O are slightly lower at the final location in this study, KoB, probably due to some local groundwater input along the Korana River course. There is no significant difference between the isotope composition of the lake water and the corresponding surface water.  Table 4, Figure 15). The SWL line for 2011-2014 is slightly steeper that the SWL for 2003-2005 (Table 4), which could be explained by the different locations and differences in climatological conditions (2003 was the warmest year in Zagreb and at the Plitvice Lakes in those periods [45]). Additionally, there is no significant difference observed between the SWL and LWL (2011-2014), although both have lower slopes than the LMWLPWLSR obtained for the PL area ( Figure 15). The lower slopes and intercepts of both SWL and LWL relations compared to the LMWL indicate the influence of the evaporation of surface waters, which is more pronounced in big lakes.   Table 4, Figure 15). The SWL line for 2011-2014 is slightly steeper that the SWL for 2003-2005 (Table 4), which could be explained by the different locations and differences in climatological conditions (2003 was the warmest year in Zagreb and at the Plitvice Lakes in those periods [45]). Additionally, there is no significant difference observed between the SWL and LWL (2011-2014), although both have lower slopes than the LMWL PWLSR obtained for the PL area ( Figure 15). The lower slopes and intercepts of both SWL and LWL relations compared to the LMWL indicate the influence of the evaporation of surface waters, which is more pronounced in big lakes.
To justify this conclusion, for all the surface and lake waters from the 2011-2014 sampling campaign we prepared the relations δ 2 H vs. δ 18 O, and the slopes and intercepts are shown in Table S4. A general decrease in both the slopes and intercepts in the downstream direction is observed. The lowest values are found in the biggest Lake Kozjak (IRB1, IRB2, and KzB).
The tritium activity concentration was determined in the surface water at the locations Matica and Lake Kozjak-Bridges in 2015, at the same time as in the groundwaters at the CR and BR springs (Figure 8, insert). The lowest mean A value was observed at CR  To justify this conclusion, for all the surface and lake waters from the 2011-2014 sampling campaign we prepared the relations δ 2 H vs. δ 18 O, and the slopes and intercepts are shown in Table  S4. A general decrease in both the slopes and intercepts in the downstream direction is observed.
The lowest values are found in the biggest Lake Kozjak (IRB1, IRB2, and KzB).
The tritium activity concentration was determined in the surface water at the locations Matica and Lake Kozjak-Bridges in 2015, at the same time as in the groundwaters at the CR and BR springs (Figure 8, insert). The lowest mean A value was observed at CR (4.3 ± 0.4 TU, n = 5), and the highest at BR (5.5 ± 0.1 TU, n = 4), and the A values the surface waters were in between them (Ma: 4.8 ± 0.2 TU, n = 5, KzB: 5.1 ± 0.2, n = 4) ( Figure S3), as one would expect after the merging of the Bijela Rijeka River and the Crna Rijeka River into the Matica River.
Although the range of δ 18 O values in the lake water is not large, systematic data can help in identifying extreme hydrological events, similar to the event in 1983-1984 explained earlier. The relation between δ 18 O and water temperature ( Figure 16) shows a generally increasing trend, as would be expected from the general relation of δ 18 O in precipitation with temperature [23]. The influence of heavy summer rains and snow melting was observed by a slight increase and decrease, respectively, in the δ 18 O values compared to the average ("normal") values in both the surface and lake waters. The higher δ 18 O data points (group a, Figure 16) are the consequences of heavy summer rains. For example, in May and June 2012 the monthly precipitation amounts of 90 mm and 144 mm were recorded in Zagreb, with δ 18 O values of −4.63‰ and −5.32‰, respectively [23]. Unfortunately, data neither for the precipitation amount nor the δ 18 O in precipitation are available for these months at the PL area. Nevertheless, the effect of the isotopically heavier precipitation is visible in group a (Figure 16), as well as in Figure S2. The isotope data of group b (Figure 16) are a consequence of the abundant isotopically lighter precipitation (rain and snow) during January-March 2013; the amount of precipitation was 148, 105, and 126 mm in January, February, and March 2013, respectively, with the δ 18 O values of −12.04‰, −13.73‰, and −11.22‰, respectively (all the data are for Zagreb [23]). The effect of the relatively lower δ 18 O and δ 2 H values was observed in all the surface and lake waters in April 2013 (Figure 14a and Figure S2).  Although the range of δ 18 O values in the lake water is not large, systematic data can help in identifying extreme hydrological events, similar to the event in 1983-1984 explained earlier. The relation between δ 18 O and water temperature ( Figure 16) shows a generally increasing trend, as would be expected from the general relation of δ 18 O in precipitation with temperature [23]. The influence of heavy summer rains and snow melting was observed by a slight increase and decrease, respectively, in the δ 18 O values compared to the average ("normal") values in both the surface and lake waters. The higher δ 18 O data points (group a, Figure 16) are the consequences of heavy summer rains. For example, in May and June 2012 the monthly precipitation amounts of 90 mm and 144 mm were recorded in Zagreb, with δ 18 O values of −4.63% and −5.32% , respectively [23]. Unfortunately, data neither for the precipitation amount nor the δ 18 O in precipitation are available for these months at the PL area. Nevertheless, the effect of the isotopically heavier precipitation is visible in group a (Figure 16), as well as in Figure S2. The isotope data of group b ( Figure 16) are a consequence of the abundant isotopically lighter precipitation (rain and snow) during January-March 2013; the amount of precipitation was 148, 105, and 126 mm in January, February, and March 2013, respectively, with the δ 18 O values of −12.04% , −13.73% , and −11.22% , respectively (all the data are for Zagreb [23]). The effect of the relatively lower δ 18 O and δ 2 H values was observed in all the surface and lake waters in April 2013 (Figure 14a and Figure S2).

Conclusions
In this paper, we presented an overview of the application of tritium and stable isotopes of water ( 2 H, 18 O) on the characterization of different water bodies (precipitation, groundwater, surface water, lake water) of the karst area of the Plitvice Lakes, Croatia. Various studies were performed over a relatively long time period (1979-2018, with a gap between 1990 and 2001). The available climatological data (amount of precipitation, air temperature) were also analyzed in a search for the evidence of climatological changes. The isotope data were compared to the continuous long-term data record for Zagreb precipitation.
The main conclusion/results of this overview are the following: • An increase in the mean annual air temperatures of 0.06 °C per year is observed for the period 1986-2019. An increase in the annual precipitation amount is also observed, at about 10 mm per year, and the range of the monthly precipitation amounts is higher in recent years.

•
Tritium activity concentration in the Plitvice Lakes precipitation shows characteristics typical for the northern hemisphere. A good correlation of the tritium activity concentration in the precipitation at the Plitvice Lakes with that in the Zagreb precipitation is observed, implying that the tritium data for Zagreb can be safely used for the study of the PL area if/when data do not exist for the PL precipitation.

•
The range of the δ 18  The deuterium excess for precipitation at the Plitvice Lakes is higher than that for the Zagreb precipitation. This is caused by the combination of higher altitude and the more intense influence of the Mediterranean air masses. The seasonal pattern in d-excess is similar to that in Zagreb, with higher values in the autumn precipitation due to the higher influence of precipitation of Mediterranean origin.

Conclusions
In this paper, we presented an overview of the application of tritium and stable isotopes of water ( 2 H, 18 O) on the characterization of different water bodies (precipitation, groundwater, surface water, lake water) of the karst area of the Plitvice Lakes, Croatia. Various studies were performed over a relatively long time period (1979-2018, with a gap between 1990 and 2001). The available climatological data (amount of precipitation, air temperature) were also analyzed in a search for the evidence of climatological changes. The isotope data were compared to the continuous long-term data record for Zagreb precipitation.
The main conclusion/results of this overview are the following: • An increase in the mean annual air temperatures of 0.06 • C per year is observed for the period 1986-2019. An increase in the annual precipitation amount is also observed, at about 10 mm per year, and the range of the monthly precipitation amounts is higher in recent years.

•
Tritium activity concentration in the Plitvice Lakes precipitation shows characteristics typical for the northern hemisphere. A good correlation of the tritium activity concentration in the precipitation at the Plitvice Lakes with that in the Zagreb precipitation is observed, implying that the tritium data for Zagreb can be safely used for the study of the PL area if/when data do not exist for the PL precipitation.

•
The range of the δ 18 O in precipitation was about 14% , and in δ 2 H it was about 110% . Various regression approaches for the determination of LMWL were applied, and the best results were obtained by the RMA and PWLSR approaches, which gave also the best results for the stations Zagreb and Zadar. The LMWL PWLSR for PL is δ 2 H = (7.97 ± 0.12) δ 18 O + (13.8 ± 1.3), n = 36, for the period 2003-2006. This LMWL lies in between the LMWL for Zagreb and Zadar.

•
The deuterium excess for precipitation at the Plitvice Lakes is higher than that for the Zagreb precipitation. This is caused by the combination of higher altitude and the more intense influence of the Mediterranean air masses. The seasonal pattern in d-excess is similar to that in Zagreb, with higher values in the autumn precipitation due to the higher influence of precipitation of Mediterranean origin.

•
The ranges in the δ values in groundwaters (at springs) were much narrower than those in precipitation, at less than 1% in δ 18 O and less than 8% in δ 2 H, indicating the good mixing of waters in karst aquifers.

•
The higher mean δ values in all three karst springs were observed in recent decades (2003-2019) than in the older period [1979][1980][1981][1982][1983][1984][1985][1986][1987][1988][1989][1990]. It can be attributed to the increase in the mean air temperature. • Different values of mean residence time (MRT) were obtained for three main springs (between 2 years for Crna Rijeka spring and 4 years for Bijela Rijeka spring). Extreme climatological and hydrological conditions in 1983-1984 enabled the estimation of the proportions of precipitation in the spring water: the shorter the MRT, the higher proportions of precipitation.

•
The amplitude in the δ 2 H of surface and lake water is less than 10% , and < 1.4% in δ 18 O, similar to that in the groundwaters. A slight increase in both the mean δ 2 H and δ 18 O values and their seasonal variations is observed for locations along the water course due to the evaporation of surface waters. No difference between lake and surface waters is observed if the mean values are concerned. There is no significant difference observed between the Surface Water Line (SWL) and Lake Water Line (LWL) (2011-2014), both having lower slopes than the LMWL PWLSR obtained for the PL area. The stable isotope composition of the surface and lake waters reacts to the extreme hydrological conditions.
The long-term and comprehensive isotope study of different water bodies in the area of the Plitvice Lakes can be an example of how the application of water isotopes ( 2 H, 18 O, 3 H) can help in the characterization of karst aquifers on the regional and global scales. The presented data have also shown what could be the topic of future research. A systematic long-term monitoring of water isotopes should be established, which could, together with the new data from the hydrogeological research, result in more detailed definition on the groundwater flow through the research area and possibly differentiate the aquifers from which the Plitvice Lakes receive the water. The usefulness of the tritium activity concentration in hydrogeological research has become recently of less importance due to the relatively constant mean value of A in precipitation during the last three decades. However, the monitoring of A in precipitation at the Plitvice Lakes may help in distinguishing global from local influences on tritium in precipitation. The further monitoring of the stable isotope data in groundwaters and surface water in relation to the global and local climate changes can give valuable data for the protection of waters, as well as tufa and biota in the area. The development of the hydrogeological conceptual model is also foreseen.  Table S1: Number of isotope data for precipitation, groundwater, surface and lake waters from the Plitvice Lakes area; Table S2: Monthly precipitation amount P, mean monthly temperature T and isotope composition of precipitation (δ 18 O, δ 2 H, deuterium excess, tritium activity concentration A) at the Plitvice Lakes; Table S3: Basic statistics (mean values and standard deviations, minimum, median and maximum values) for δ 18 O (shadowed rows) and δ 2 H in lake waters (IRB1, IRB2, IRB3, IRB4) and surface waters, 2011-2014 period; Table S4