*Article* **Long-Term Isotope Records of Precipitation in Zagreb, Croatia**

**Ines Krajcar Broni´c 1,\*, Jadranka Bareši´c 1, Damir Borkovi´c 1, Andreja Sironi´c 1, Ivanka Lovrenˇci´c Mikeli´c <sup>1</sup> and Polona Vreˇca <sup>2</sup>**


Received: 29 November 2019; Accepted: 10 January 2020; Published: 14 January 2020

**Abstract:** The isotope composition of precipitation has been monitored in monthly precipitation at Zagreb, Croatia, since 1976. Here, we present a statistical analysis of available long-term isotope data (3H activity concentration, δ2H, δ18O, and deuterium excess) and compare them to basic meteorological data. The aim was to see whether isotope composition reflected observed climate changes in Zagreb: a significant increase in the annual air temperature and larger variations in the precipitation amount. Annual mean δ18O and δ2H values showed an increase of 0.017% and 0.14% per year, respectively, with larger differences in monthly mean values in the first half of the year than in the second half. Mean annual *d*-excess remained constant over the whole long-term period, with a tendency for monthly mean *d*-excess values to decrease in the first half of the year and increase in the second half due to the influence of air masses originating from the eastern Mediterranean. Changes in the stable isotope composition of precipitation thus resembled changes in the temperature, the circulation pattern of air masses, and the precipitation regime. A local meteoric water line was obtained using different regression methods, which did not result in significant differences between nonweighted and precipitation-weighted slope and intercept values. Deviations from the Global Meteoric Water Line GMWL (lower slopes and intercepts) were observed in two recent periods and could be explained by changes in climate parameters. The temperature gradient of δ18O was 0.33%/◦C. The tritium activity concentrations in precipitation showed slight decreases during the last two decades, and the mean *A* in the most recent period, 2012–2018, was 7.6 ± 0.8 Tritium Units (TU).

**Keywords:** precipitation; Zagreb; Croatia; stable isotope ratios; 2H/ 1H and 18O/ 16O; deuterium excess; local meteoric water line; δ18O–temperature relation; tritium

#### **1. Introduction**

Water, especially groundwater, has become an invaluable natural resource, and the availability of freshwater is one of the greatest issues facing mankind today [1]. A consistent and careful assessment and management of water resources is crucial for their sustainable development. This can be performed by various methodologies, among which isotope methods using environmental (stable and radioactive) and artificial radioactive isotopes have proven to be effective tools for solving many critical hydrological problems and processes [2–9]. In many cases, isotope techniques have provided information that could not be obtained through any other conventional means [10–12].

Isotopes that are constituent elements of a water molecule are of special interest as perfect candidates for water tracers: hydrogen (1H, 2H, 3H) and oxygen (16O, 17O, 18O) isotopes. Among these, only 3H is a radioactive isotope, while the others are stable isotopes.

Precipitation presents an input to groundwater, and therefore knowledge on the isotope composition of precipitation is a prerequisite for groundwater studies. Temporal and spatial patterns of isotopes in precipitation (expressed as δ2H and δ18O values and the tritium activity concentration, *A*) have been observed since the 1950s and have contributed to hydrological and hydrogeological research [13–16], climate and paleoclimate studies [4,15,17–22], and ecological research [23–25]; further, precipitation isotope mapping has been widely implemented during recent decades [25,26].

The International Atomic Energy Agency (IAEA) and the World Meteorological Organization (WMO) have recognized the importance of the isotopic composition of precipitation on a global scale. A program involving the worldwide monitoring of the isotopic composition of monthly precipitation (called the Global Network of Isotopes in Precipitation (GNIP)) was therefore established in 1961. The objective of the network was a systematic collection of data on the isotopic composition of precipitation across the globe to determine temporal and spatial variations of isotopes in precipitation. Isotopic data include the tritium activity concentration, *A* (expressed in tritium units, TU), and stable isotopes of hydrogen and oxygen isotopes (δ2H and δ18O values), as well as climatological data (mean monthly temperature, monthly precipitation amount, and atmospheric water vapor pressure) [27]. The collected records have enabled the establishment of seasonal variations and various correlations among the data. Seasonal variations in δ18O and δ2H values of precipitation and their weighted mean annual values have remained fairly constant from year to year at a given location as long as the annual range and sequence of climatic conditions did not change significantly from year to year.

However, in recent years, we have become aware of climatic changes that have caused an increase in global temperature and changes in the precipitation pattern, as well as severe and extreme weather events (droughts, heavy storms, and temperature records). The recent increase in global temperature has exceeded the natural variabilities during the Holocene [28]. The planet's average surface temperature has risen about 1 ◦C (between 0.8 ◦C and 1.2 ◦C [29]) since the late 19th century, a change driven largely by increased carbon dioxide and other human-made emissions into the atmosphere. Most of the warming has occurred in the past 35 years, with 16 of the 17 warmest years on record occurring since 2001. The 10 hottest years ever recorded have all occurred since 1998 [30]. The hottest years on record globally have been the last five (2014–2018), with 2016 being the hottest year [31], and eight months in 2016 (from January through September, except for June) were the warmest on record for those respective months. October, November, and December of 2016 were the second warmest of those months on record: in all three cases, behind records set in 2015. Correlations between the precipitation isotope ratios recorded in the GNIP and meteorological quantities may provide additional evidence of recent climate change that appears to have manifested globally as well as evidence of the local weather situation. To find such evidence, one should have sufficiently long records of both climate data and the isotopic composition of precipitation.

Monitoring of the isotope composition of monthly precipitation at a station in Zagreb (Croatia) has been performed since 1976 (tritium activity concentration, *A*) and since 1980 (stable isotope ratios of hydrogen (2H/ 1H) and oxygen (18O/ 16O)). Isotope data up to 2003 are available in the GNIP database [27]. This work presents details on the history of monitoring the isotope composition of precipitation in Zagreb, Croatia, for the period 1976–2018. Such series of isotope data are rather scarce in Europe [27], and the present time series analysis can be a first step toward a more detailed comparison of data reported for sites with similar long-term records, such as Vienna (Austria), Krakow (Poland), and Ljubljana (Slovenia), in order to obtain a wider spatial and temporal pattern of the isotope composition of precipitation. Statistical analyses of isotope (δ2H, δ18O, deuterium excess, and tritium activity concentration) and basic meteorological data (temperature and precipitation amount) were performed. The complete 43-year-long record was divided into subperiods in order to better investigate climatological and isotope-in-precipitation changes. Subperiods should be neither too long nor too short, so we arbitrarily chose four almost equally long subperiods: 1976–1985, 1986–1995, 1996–2006, and 2007–2018.

In the following section, we introduce the notation (δ2H, δ18O) used for stable isotope composition as well as the concepts of meteoric water lines and deuterium excess. We describe the behavior of tritium activity concentration in the atmosphere and give a brief description of measurement techniques that have changed during the studied period. An overview of sampling locations of monthly precipitation in Croatia and climate characteristics of the area are presented, which will help in discussing the data from Zagreb. In Section 3, we show the results of the monthly data, while in Section 4, we present a discussion of the statistical analyses of the annual data and average values in the subperiods, the observed temporal trends, and various correlations among the data. We discuss how the observed temperature and precipitation amount changes were recorded in the isotopic composition of precipitation in Zagreb.

#### **2. Materials and Methods**

#### *2.1. Stable Isotopes of Hydrogen and Oxygen*

The isotopic composition of water constituents depends on isotope fractionation caused by phase transfers of water masses (evaporation/precipitation), which depend on the area of water origin (latitude, altitude, continent or maritime, climate region) and the precipitation amount [2,4,32–34]. Therefore, the isotopic composition of precipitation of different origins and seasons provides for the application of stable water isotopes as tracers of the hydrological cycle.

The results are reported as δ-values per mill (%) relative to the standard [4,35–38]:

$$
\delta\_{\rm S/R} = \frac{R\_{\rm Sample}}{R\_{\rmReference}} - 1.\tag{1}
$$

Here, *R*Sample and *R*Reference stand for the isotope ratio (*R* = 2H/ 1H and *R* = 18O/ 16O) in the sample and the reference material (standard), respectively. Standard mean ocean water, SMOW, has been proposed as a (virtual) standard for reporting measured values [39]. SMOW is an arbitrary mean value based on the Epstein–Mayeda oxygen scale obtained from deep ocean water, since it does not interact with the atmosphere and has a stable isotopic composition [40], and it is defined in terms of an actual water reference standard, the NBS-1 (National Bureau of Standards, USA). In 1968, the IAEA established an international standard, the Vienna SMOW (VSMOW), which has been replaced by the VSMOW2 [38,41].

The δ2H and δ18O isotopic compositions of meteoric waters (precipitation and atmospheric water vapor) are strongly correlated. If δ2H is plotted versus δ18O, the data cluster along a straight line is

$$
\delta \, \delta^2 \mathcal{H} = 8.0 \cdot \delta^{18} \mathcal{O} + 10. \tag{2}
$$

The relation in Equation (2) is referred to as the global meteoric water line (GMWL) [2,39,42]. It describes the general relation between δ2H and δ18O on a global scale reasonably well. However, the intercept is higher for precipitation originating from the Mediterranean area [2,4,43–45], while the slope does not change, as in the case of the eastern Mediterranean meteoric water line, δ2H = 8 δ18O + 22 [43]. This example shows that for applications in hydrogeological studies, regional local meteoric water lines (LMWLs), either long-term or for certain shorter periods, can be more appropriate. Generally, an LMWL has the form δ2H = *a* δ18O + *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 [1,4,46,47]. There are different ways of calculating *a*- and *b*-values. Traditionally, the ordinary least squares regression (OLSR) is used, and more recently, the reduced major axis regression (RMA) and the major axis least squares regression (MA, sometimes also called the orthogonal regression) [48,49] have been applied. The precipitation-weighted OLSR approach (PWLSR) was introduced to reduce the impact of small precipitation events that are hydrogeologically not

significant [49]. Similarly, precipitation-weighted RMA and MA, i.e., PWRMA and PWMA, regressions have been applied to data from the GNIP database, with at least 36 monthly datapoints available [48].

Deuterium excess (*d*-excess, or *d*) is defined as [2]

$$d = \delta^2 \mathcal{H} - 8 \,\delta^{18} \mathcal{O}\_{\prime} \tag{3}$$

which can be related to the meteorological conditions in the source region from which the water vapor is obtained [1,4,35,43,46]; therefore, it can be used to identify vapor source regions, and it is often considered to be the most useful parameter in characterizing vapor origin [43]. Winter precipitation originating from the Mediterranean Sea is characterized by distinctly higher *d*-excess values (*d* > 18%) than is precipitation coming from the Atlantic (*d* ~ 10%), reflecting the specific source conditions during water vapor formation. Increased deuterium excess in precipitation can also arise from a significant addition of re-evaporated moisture from continental basins to water vapor traveling inland [1,4,27,35,44–46,50,51].

#### *2.2. Tritium*

Tritium (3H) is a natural cosmogenic isotope of hydrogen that is formed in the upper atmosphere through reactions of thermal neutrons with 14N. It oxidizes to tritiated water, H3HO, and thus enters the natural water cycle. The half-life of 3H is 12.32 years [52], and it decays to 3He by emitting beta particles with a maximal energy of 18.6 keV.

Tritium is also an anthropogenically produced isotope, and it can be differentiated as "bomb-produced" tritium and technogenic tritium. Massive injections of 3H from weapons tests in the 1950s and 1960s, mostly in the Northern Hemisphere, caused an almost 100-fold increase in the tritium activity concentration in precipitation [1], known as the bomb peak. The highest concentration of tritium, about 6000 TU (1 TU = 0.118 Bq/l), was observed in 1963 in precipitation at continental stations in the Northern Hemisphere [1,27,53], while at maritime stations, the maximal values were lower (about 2000 TU). The data at the marine stations were systematically lower than at the continental stations because moisture evaporated from the ocean has a low 3H activity concentration due to the long residence time of water in the ocean. After the cessation of atmospheric nuclear weapons tests, a gradual decrease in 3H activity concentration at all stations in both hemispheres was observed due to natural decay and the washout of tritium into the oceans and groundwater. The levels of tritium have declined globally and regionally, approaching the natural pre-bomb level. The pre-bomb natural tritium activity concentration is assumed to be about 1 TU in oceanic regions, about 10 TU in inland areas, about 5 TU in central Europe [1], and about 5 TU on average globally [54]. Monitoring of the tritium level in precipitation at several short-distance stations showed that there was no significant systematic discrepancy between them [1]. The "anthropogenically modified natural distributions" present now are "new natural global" environmental levels.

Technogenic tritium is produced in various industries, such as nuclear power plants, nuclear reactors, future fusion reactors, fuel reprocessing plants, heavy water production facilities, medical diagnostics, radiopharmaceuticals, luminous paints, sign illumination, self-luminous aircraft, airport runway lights, luminous dials, and gauges and wrist watches [55–58]. Technogenic tritium causes deviations from the "anthropogenically modified natural tritium distribution" at a local or regional level.

The seasonal and spatial distribution of tritium activity concentration in precipitation around the globe has been found to be dominated by the annual stratosphere–troposphere exchange at high latitudes in early spring, in combination with latitudinal and continental effects [1]. The latitude effect is described as the highest 3H activity concentrations observed between the 30th and 60th parallel, with values lower by a factor of approximately five at low-latitude and tropical stations.

It should be noted that the seasonal variations in 3H activity concentration in precipitation do not have the same origin as the seasonal variations in δ2H and δ18O: variations in 3H are caused by the exchange between the stratosphere and the troposphere and are not caused or influenced by the

temperature, while the local seasonal variations in δ2H and δ18O show a close relation with the local temperature [1,2,4,47].

#### *2.3. Sampling Sites and Climate*

There are three main climate types prevailing in Croatia: continental, maritime, and mountain. Such a climate distribution is determined by the geographical position of Croatia in the northern midlatitudes and the corresponding weather processes. Croatia is a relatively small country (56,594 km2) positioned between the Pannonian Plain and the Adriatic Sea and has a large orographic variety. Therefore, the most important climate modifiers are the Adriatic Sea in the southwest, the mountain chain Dinarides in the central part, and openness to the Pannonian plain in the northeast [59]. Accordingly, most of Croatia has a temperate rainy climate (Köppen code Cxx [60,61]). For example, the Zagreb climate zone is described as Cfb: a temperate climate without a dry season and with a warm summer. The highest mountain areas have a cold snow and forest climate (Köppen code Dxx). The complete designation of climates for particular sampling sites is presented in Table S1 (Supplementary Materials). A comparison of the 30-year period 1981–2010 to the standard climatological period, 1961–1990, showed a significant increase in the mean annual temperature at all 20 studied stations in Croatia, with stronger warming at the continental stations than along the coast and with the largest changes in the summer [59]. It was also noted that both minimal and maximal temperatures increased by larger amplitudes inland compared to along the coast. The precipitation pattern has also changed, but differently in different parts of Croatia: an increase in the precipitation amount has been observed inland, with a statistically significant increase in autumn [59]. The observed climate changes (temperature increase and the precipitation regime) have resulted in changes in climate classes for some stations. The stations with isotope-in-precipitation data for which the climate class changed are Dubrovnik and Zadar (from Cfa to Csa) and Puntijarka on Mt. Medvednica, near Zagreb (from Dfb to Cfb) (Table S1).

A long-term isotopes-in-precipitation record (1976–2018) exists in Croatia only for the Zagreb station [27,47,62] (although microlocations have changed, as will be explained later). Data for some other stations with shorter monitoring periods were obtained during various individual projects (location numbers 1 to 12, Figure 1) [18,47,62–68]. Details on sampling sites are presented in Table S1. Some projects have also included monitoring at stations in Slovenia (location numbers 13 to 15, Figure 1), e.g., in Ljubljana [47,62,69–71], Portorož, and Kozina [47,72,73].

Precipitation monthly composite samples for the Zagreb station were collected at the Ruđer Boškovi´c Institute (RBI, 45.817◦ N, 15.967◦, 165 m a.s.l.) from 1976 until 1995. In 1994 and 1995, higher 3H activity concentrations were measured due to experimental research in the nearby department in which technogenic tritium was used [74,75], while there was no influence on stable isotope data. As a consequence of local tritium contamination, a new location for precipitation sampling had to be found. During 1995 and 1996, precipitation samples were additionally collected at Puntijarka, on Mt. Medvednica (45.917◦ N, 15.95◦, 988 m a.s.l.), 15 km north of the city of Zagreb [62], and at the Zagreb–Griˇc site at the Croatian Meteorological and Hydrological Service (45.814◦ N, 15.972◦ E, 157 m a.s.l.) in the center of Zagreb (in 1996). 3H activity concentrations in precipitation at the stations Zagreb–Griˇc and Puntijarka in 1996 were almost identical, 10.7 ± 4.0 TU and 10.4 ± 6.0 TU, respectively [62]. For the analyses in this paper, the Zagreb–RBI data for tritium activity concentration were used for up to 1993, data from Puntijarka were taken as representative for Zagreb in 1995, and for 1996 and on, data from Zagreb–Griˇc were used. Tritium data exist for 1994, but they will not be discussed here because they present technogenic (local) tritium. Data for the stable isotope composition of Zagreb precipitation are not available for the 2007–2009 period and for 2011.

Daily rain events were collected in Zagreb for the period from October 2002 to March 2003, and the obtained stable isotope composition was compared to the monthly δ18O and δ2H data [76].

**Figure 1.** Map of stations with isotope-in-precipitation data in Croatia (locations 1–12) and in Slovenia (locations 13–15). Details on stations are in Table S1.

#### *2.4. Meteorological Data*

The meteorological data consisted of monthly precipitation amount and average monthly air temperature. Data were obtained on request from the Croatian Meteorological and Hydrological Service (CMHS). Meteorological records for Zagreb, Croatia, exist for the period since 1862 and have been analyzed by CMHS for the usual climatological periods [77]. Here, we used data from only the 1976–2018 period, for which we had records of the isotope composition of precipitation. Minimal and maximal monthly values within a year were identified, and the mean annual values of temperature and the total annual precipitation amount were determined.

#### *2.5. Measurement of* δ*2H and* δ*18O*

The stable isotope composition of the precipitation samples for the period 1980–2003 was measured on a Varian MAT 250 dual inlet isotope ratio mass spectrometer (IRMS) at the Jožef Stefan Institute in Ljubljana [49,69,70,73]. The isotopic composition of hydrogen (δ2H) was determined by means of the H2 generated by the reduction of water over hot zinc (up to 1998) [69], and later over hot chromium [78]. The oxygen isotopic composition (δ18O) was measured by means of the water–CO2 equilibration technique [40]. All measurements were carried out together with laboratory standards that were calibrated periodically against international standards, as recommended by the IAEA. The measurement precision of duplicates was better than <sup>±</sup>0.1% for <sup>δ</sup>18O and <sup>±</sup>1% for <sup>δ</sup>2H. The stable isotope composition of the precipitation sampled between 2004 and 2006 was determined at SILab (Stable Isotope Laboratory at the Physics Department, School of Medicine, University of Rijeka, Rijeka, Croatia). An HDO Equilibration Unit (ISO Cal) attached to the dual inlet port of a DeltaPlusXP (Thermo Finnigan) IRMS was used [68,79]. The δ18O and the δ2H were obtained from CO2 and H2 gas, respectively, after equilibration with a 4-ml water sample. The measurement reproducibility of duplicates was better than <sup>±</sup>0.1% for <sup>δ</sup>18O and <sup>±</sup>1% for <sup>δ</sup>2H. The stable isotope composition

of precipitation from 2010 and on was analyzed with a Liquid Water Isotope Analyzer (LWIA-24d, Los Gatos Research) at the Institute for Geochemical Research, the Hungarian Academy of Sciences, Budapest, Hungary. The uncertainty of the measurements was reported to be <sup>±</sup>0.2% for <sup>δ</sup>18O and <sup>±</sup>0.6% for <sup>δ</sup>2H [66]. The stable isotope composition of Zagreb's precipitation for the 2012–2018 period 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). Data were analyzed by the Laboratory Information Management System (LIMS) [80]. The measurement precision of duplicates was <sup>±</sup>0.1% for <sup>δ</sup>18O and <sup>±</sup>0.3% for <sup>δ</sup>2H.

#### *2.6. Measurement of Tritium Activity Concentration*

Tritium activity concentration (*A*) in all monthly samples was determined at the Ruđer Boškovi´c Institute in Zagreb. The results are expressed in tritium units (1 TU = 0.118 Bq l<sup>−</sup>1) [1] since the same units are used by the IAEA–WMO/GNIP database [27]. A tritium unit represents one 3H atom in 10<sup>18</sup> atoms of hydrogen. The gas proportional counting technique (GPC) was used up to 2009 [81–83]. Methane (CH4) was used as a counting gas in a multiwire proportional counter. It was obtained through the reaction of water (50 ml) with aluminum carbide at 150 ◦C [81]. The counting energy window was set to energies between 1 keV and 10 keV to obtain the best figure of merit. Gas quality control was performed by simultaneously monitoring the count rate above the tritium channel, i.e., above 20 keV [82]. The detection limit was 2.5 TU, and the measurement uncertainty was between 2 and 5 TU, depending on the activity concentration. In 2008, a technique of liquid scintillation counting of electrolytically enriched samples (LSC-EE) was introduced, and between 2008 and 2009, GPC and LSC-EE techniques were used [83,84]. Since 2010, samples have been measured using the LSC-EE technique only [71,83–86].

The electrolytic enrichment system at the RBI was produced by the AGH University of Science and Technology, Krakow, Poland [87,88]. It consists of 20 cells 500 ml in volume (stainless steel anodes and mild steel cathodes). Each sample was distilled before electrolysis, as the required conductivity is <50 μS/cm. An enrichment run contained 15 unknown samples, 3 spike waters (water of known tritium activity concentration, 500–600 TU) used for monitoring the electrolysis performances and the calculation of enrichment factor, and 2 tritium-free samples used for system control. The enrichment procedure at RBI took 1420 Ah distributed over 8 days. Enriched samples were distilled again after electrolysis, with 6–8 g of PbCl2 added to each sample. The scintillation cocktails for measurement in LSC were prepared by 8 ml of sample and 12 ml of scintillator Ultima Gold LLT in high-density low-diffusion polyethylene vials. Measurements were performed by an ultra-low-level liquid scintillation counter, the Quantulus 1220, in 10 cycles of 50 min each. Each measurement run consisted of 24 scintillation cocktails: 20 enriched samples, 1 nonenriched spike sample, 1 international standard, and 2 nonenriched background samples. On the basis of the initial and final mass of water in cells and the individual count rates of the spike water before and after enrichment, the enrichment factor *E* was calculated [86,87]. The average 10-year *E* value of the system was 26 ± 2. The detection limit obtained by the LSC-EE technique was around 0.5 TU, and the measurement uncertainty was between 0.5 and 3 TU [71,83–86].

#### *2.7. Data Evaluation*

Mean δ18O, δ2H, and *d*-excess values, weighted by precipitation amount, were calculated from all monthly data and then summed over all collected samples per month and per year. Thus, the monthly mean values for a specific month over a certain period were obtained, as were the annual mean values. The number of datapoints per year for all years (Table 1) fulfilled the requirement for the calculation of annual mean values, i.e., the lowest number of monthly isotope datapoints per year was 9, while at least 7 monthly samples were required [89]. All years with an incomplete number of monthly samples, however, met the requirement that the available isotope data comprise more than 70% of the total precipitation amount collected per year (Table 1) [89].


**Table 1.** Mean annual temperature *T*, ranges of monthly temperatures within a year, annual precipitation amount *P,* weighted mean (w.m.) annual δ18O, δ2H, *d*-excess (*d*), and the mean annual tritium activity concentration *A* in precipitation at Zagreb, 1976–2018. Here, *n*: number of monthly datapoints; %*P*: percentage of precipitation if different from 100%, comprised by δ18O, δ2H, and *d*-excess data; TU: tritium unit. Bold font: the lowest and highest values in a series.

Correlations between various datapoints were obtained as ordinary least squares regressions, and the Pearson's coefficient *r* is given, as are the number of data pairs *n* and the *p*-value describing the statistical significance of the correlations. Data taken from the literature usually have the adjacent *r*<sup>2</sup> value reported.

In the special case of correlations between δ2H and δ18O (i.e., for LMWLs), different methods were applied: an ordinary least squares regression (OLSR), a reduced major axis regression (RMA), and a major axis least squares regression (MA) [48,49,89]. In addition, we calculated precipitation-weighted regressions (PWLSR, PWRMA, and PWMA) [48,49], which took into account the precipitation amount in a particular month. The local meteoric water lines are defined as LMWLOLSF, LMWLRMA, LMWLMA, LMWLPWLSR, LMWLPWRMA, and LMWLPWMA. While OLSR regressions were found by commercially available software (MS Excel), for other regressions we used Local Meteoric Water Line Freeware [90]. 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 for that set of data [48].

Deuterium excess *d* was calculated from paired monthly data according to Equation (3). Some precipitation samples showed very low (highly negative; in our case, the lowest value was −13%) *d*-excess values. This can be caused by improper sampling (e.g., the precipitation stayed for the whole month in a sample collector), by evaporation due to a low amount of precipitation and/or high temperatures, or by the evaporation/sublimation of raindrops falling in a dry atmosphere [50,51,71].

#### **3. Results**

Data on monthly temperatures, monthly precipitation amount, and the tritium activity concentration of precipitation at Zagreb for the 1976–2018 period, as well as the stable isotope composition of monthly precipitation (δ18O, δ2H, *d*-excess) for the 1980–2018 period, are shown in Table S2 (Supplementary Materials).

#### *3.1. Meteorological Data*

Minimal monthly (mean air) temperatures at Zagreb for the 1976–2018 period ranged from −3.4 ◦C (in 1985) to 5.4 ◦C (in 2014) (Figure 2a). They were measured in January (in 21 cases, or 49%), February (12 cases, or 27.9%), and December (9 cases, or 20.9%), and only once, in 1988, was the coldest month November. The highest monthly temperatures, in a range from 18.9 ◦C (1978) to 25.8 ◦C (1992, 2003) were measured in 26 out of 43 years (60.5%) in July, 15 times (35%) in August, and only in 2 cases in June (1979 and 1996).

**Figure 2.** *Cont*.

**Figure 2.** (**a**) Mean monthly air temperatures; (**b**) monthly precipitation amount. Data from station Zagreb–Griˇc for the 1976–2018 period. Data obtained from the Croatian Meteorological and Hydrological Service (CMHS).

Mean annual temperatures, together with yearly minimal and maximal monthly temperatures and the total yearly amount of precipitation at Zagreb for the 1976–2018 period, are shown in Table 1. Mean annual temperatures ranged from 9.6 ◦C in 1980 to 14.1 ◦C in 2018.

The annual precipitation amount ranged from 521 mm (2011) to 1234 mm (2014) (Table 1), while the maximal monthly precipitation amount occurred in August 1989 (260 mm), followed by 236 mm in September 2017 and 208 mm in September 2014. In all other months, the monthly amount of precipitation was below 200 mm (Figure 2b).

#### *3.2. Stable Isotopes*

Monthly <sup>δ</sup>18O values in precipitation at Zagreb ranged from <sup>−</sup>17.6% (January 2005) to <sup>−</sup>0.5% in June 1998 (Figure 3). Similarly, the lowest <sup>δ</sup>2H = <sup>−</sup>133.1% (January 2005) and the highest <sup>δ</sup>2H = −11.4% (June 1998) were determined in the same months (Figure S1, Supplementary Materials).

The lowest weighted mean annual values were observed in 2010 (δ18O <sup>=</sup> <sup>−</sup>9.7%, <sup>δ</sup>2H <sup>=</sup> <sup>−</sup>68.0%) and the highest in 2000 (δ18O = <sup>−</sup>5.54%, <sup>δ</sup>2H = <sup>−</sup>39.68%) (Table 1). The number of monthly isotope datapoints (*n*) available in each year as well as the percentage of annual precipitation amount (%*P*) comprised by the isotope data during the *n* months (Table 1) satisfied the requirements for the calculation of mean annual values [89]. For 15 years out of 39 (with the available stable isotope data), the number of monthly samples was less than 12, but it was never less than 9, and the available isotope data comprised at least 77% of annual precipitation in these years.

**Figure 3.** Monthly δ18O in precipitation at Zagreb, 1980–2018 period.

Monthly values of deuterium excess (Figure 4, Table S2) ranged from −6.9% in May 2000 to 22.5% in October 1994, with 110 mm rain, probably from the Mediterranean. It should be noted here that for six months of the entire studied period, the monthly *d*-excess values were lower than −7%, i.e., they were more than three standard deviations lower than the overall mean *d*-excess value and were therefore deleted from the record and further analysis. Three out of six cases were caused by the evaporation of small monthly precipitation amounts (11 mm in March 1996, 22 mm in February 2000, and 2 mm in December 2016), and the other three (January 1996, June 1998, and July 2012) were probably from the evaporation/sublimation of raindrops. Mean annual *d*-excess values ranged between 2.28% and 10.74% in 2000 and 1995, respectively (Table 1).

**Figure 4.** Monthly deuterium excess values for precipitation at Zagreb, 1980–2018.

#### *3.3. Tritium Activity Concentration in Precipitation at Zagreb*

The complete record of tritium activity concentration (*A*) in precipitation for Zagreb, 1976–2018 (Figure 5), exhibited a pattern typical of continental stations of the Northern Hemisphere. Seasonal variations were superposed on the basic decreasing trend of mean annual values until approximately 1996. The maximal monthly 3H activity concentration at the Zagreb station was observed between May and July, mostly in June. A secondary maximum was also observed three times in January and February. The lowest 3H activity concentrations were almost uniformly distributed from October to February, with a slightly more frequent occurrence in December.

**Figure 5.** Complete record of monthly tritium activity concentration (*A*) in precipitation at Zagreb, 1976–2018. Insert: for the 1995–2018 period.

Data recorded from 1995 to 2018 (insert in Figure 5) showed an almost constant mean annual 3H activity concentration ranging between 5.4 TU (in 2004) and 13.5 TU (in 1995), with a mean value of 8.5 ± 1.2 TU. Seasonal variations remained observable, with winter activities close to the natural pre-bomb 3H activity concentrations (≤5 TU) and summer values up to 21 TU [91,92].

#### **4. Discussion**

#### *4.1. Trends in Meteorological Parameters*

The annual precipitation amount *P* at Zagreb for the 1976–2018 period (Figure 6) showed a slight increase (1.4 ± 1.7 mm/y, *r* = 0.13, *p* = 0.4), as did the maximal monthly values within a year (0.7 ± 0.5 mm/y, *r* = 0.23, *p* = 0.14), while the minimal monthly values within a year, including no-rain months, showed a slight decrease (−0.1 ± 0.1 mm/y, *r* = −0.13, *p* = 0.41). However, the trends (statistically not significant) were not the most prominent characteristics of the data. Higher dispersion/fluctuations from the mean value for the whole period, 1976–2018 (867 ± 138 mm), were obvious (Figure 6). In the period 1976–2000, practically all values lay within ±1 standard deviation (±1 σ), while later on there were some years with deviations from the mean over ±2 σ. The mean values of the precipitation amount in the subperiods (1980–1985, 1986–1995, 1996–2006, and 2012–2018) for which isotope data were available showed larger fluctuations in the precipitation amount in more recent periods (Table 2).

**Figure 6.** Annual precipitation amount *P* in Zagreb, 1976–2018 period. Solid line: the mean value for the whole period; dashed lines: ±1σ; dotted lines: ±2σ.

**Table 2.** Comparison of mean values of meteorological parameters (*T*, *P*) and isotopic data (δ18O, δ2H, and *d*-excess) for different periods.


Mean annual temperature and the minimal and maximal monthly temperature within a year (Figure 7) showed a significant increase at a 95% confidence level (*p* < 0.05) (a mean value of 0.071 ± 0.008 ◦C per year (*r* = 0.82), a minimal value of 0.05 ± 0.08 ◦C/year (*r* = 0.33), and a maximal value of 0.09 ± 0.02 ◦C per year (*r* = 0.69)). Since both the minimal and maximal monthly temperatures increased, but with different gradients, the result was that the increase in the amplitudes of the air temperatures (0.04 ± 0.03 ◦C per year, *r* = 0.26, *p* = 0.09) was significant at a 90% significance level (Figure 7).

It is also interesting to look at the monthly mean *P* and *T* values at the Zagreb station averaged for each month over a certain period (Figure 8) and calculated for four subperiods. The monthly amount of precipitation was relatively uniformly distributed throughout the year. However, the average values in January–April and in December were lower (<61 mm) than in May–November (>70 mm). The month with the highest average precipitation was September, with 97 mm. September was also the only month showing a constant increase in the amount of precipitation over the studied periods. Monthly precipitation for the period 1976–1996 (i.e., the first two periods from Figure 8a) also showed lower precipitation (<60 mm) for January to April. However, maximum precipitation amounts were observed in June (94 mm) and August (91 mm) [62]. Such a shift in the precipitation regime is in accordance with observations of climate changes in Croatia [59].

**Figure 7.** Mean annual temperature, minimaland maximal monthly temperatures within a year, and the temperature amplitudes in Zagreb, 1976–2018 period. Trend lines for each dataset are shown in the same color.

Mean monthly temperatures (Figure 8b) for all months in the 1976–1985 period were lower than in the most recent period, 2007–2018. The difference ranged from 1.1 ◦C in October and December to 3.5 ◦C in August and 3.6 ◦C in April. This observation once more corroborates observations of constant recent temperature increases that are more pronounced in the spring–summer periods [59].

#### *4.2. Trends in Stable Isotope Data*

Both the arithmetic mean annual values (δ18Oa, δ2Ha, *d*a) and the weighted mean annual values by amount of precipitation (δ18Ow.m., δ2Hw.m, *d*w.m.) were calculated (Table S3). Due to the relatively homogeneous distribution of the annual precipitation amount (Figure 8a), there was no significant difference between the two types of annual means; in fact, a very good correlation was obtained, and here we present an example of the δ18O values:

$$
\delta^{18} \mathcal{O}\_{\text{w.m.}} = (1.006 \pm 0.16) \,\delta^{18} \mathcal{O}\_{\text{a}} + (0.48 \pm 1.4), n = 34, r = 0.73. \tag{4}
$$

Changes in the weighted mean annual values of δ18O, δ2H, and *d*-excess in the studied 1980–2018 period are shown in Figure 9 together with their respective trends. The δ18O and δ2H values exhibited increases, with an increase rate of 0.017% <sup>±</sup> 0.014% per year (*r* = 0.21, *p* = 0.23) for <sup>δ</sup>18O and 0.14% <sup>±</sup> 0.11% per year (*r* = 0.23, *p* = 0.19) for <sup>δ</sup>2H. The annual mean *d*-excess remained constant (slope <sup>≈</sup> 0). The corresponding weighted mean values for the periods 1980–2006 (as well as for the shorter subperiods) and 2012–2018 are shown in Table 2. Both the δ18O and δ2H values were higher in the more recent period, 2012–2018, than from 1980 to 2006. We observed a similar trend earlier: the mean <sup>δ</sup>18O in the 2001–2003 period (−8.3%) was more positive than the long-term mean <sup>δ</sup>18O (−8.8%) [76]. However, it was noticed (Table 2) that the values in the 1996–2006 and 2012–2018 periods were practically the same, although the temperature differed in these periods.

(**a**)

**Figure 8.** (**a**) Monthly mean precipitation amounts for the four subperiods and the average value for the whole period; (**b**) monthly mean air temperature for the subperiods.

(**b**)

**Figure 9.** Temporal changes of weighted mean annual δ18O, δ2H, and *d*-excess values.

To investigate which months contributed most to the changes in the mean annual values of δ18O and δ2H, we calculated the monthly mean δ18O values in the four subperiods (Figure 10). The most pronounced differences were observed in the first half of the year, with higher values in the 2012–2018 period than in earlier years in January, March, and April, while in February, the most recent δ18O and δ2H were the lowest. The differences in the second half of the year were not large. This is behavior similar to that described earlier for the monthly mean temperature (Figure 8b). However, while the summer temperatures (June, July, and August) were the highest in the most recent period, the δ18O and δ2H were not.

**Figure 10.** (**a**) Monthly mean values of δ18O in the four subperiods; (**b**) monthly mean values of δ2H in the four subperiods**.**

No significant correlation was observed between δ18O and *P*. This was expected for the midlatitude continental station of Zagreb, which had an expressed seasonality in temperature (Figure 8b) and δ18O values (Figure 10a) and a lack of clear seasonality in precipitation amount (Figure 8a) [4]. No amount effect was observed or reported for stations in Croatia and Slovenia [17,47,62,70–73].

#### *4.3. Deuterium Excess*

In previous analyses [47,76], it was demonstrated that the deuterium excess values for Zagreb's precipitation were higher in autumn than in the spring. The same occurred in this study using monthly mean values of *d*-excess for the four periods, 1980–1985, 1986–1995, 1995–2006, and 2012–2018 (Figure 11). In all periods, values from January to June (<8%) were lower than those in the second half of the year (>8%). While in the first half, there was a slight decrease in the newer periods (7.6 ± 0.6% in 1980–1985 to 6.0 ± 1.6% in 2012–2018), the *d*-excess values in the second half of the year showed an increasing trend. Such an interplay of decreases and increases of the monthly mean values eventually resulted in no change in the mean annual *d*-excess values over the whole studied period, as was shown earlier (Figure 9, Table 2). The higher *d*-excess in autumn (Figure 11) indicates a higher influence of the Mediterranean air masses in these months. It is interesting to note the shift of the autumn peak in the *d*-excess value from October (in 1980–1985) to November in the 2012–2018 period. Such behavior resembled the changes in the monthly precipitation amount in autumn (Figure 8a). November was the only month in which a significant increase in *d*-excess in the whole period was observed, with a rate/slope of 0.14% ± 0.04% per year (*n* = 35, *r* = 0.54, *p* < 0.05).

**Figure 11.** Monthly mean values of *d*-excess in the four subperiods.

A similar pattern of monthly *d*-excess distribution was also observed for all other Croatian stations (Figure 1, Table S1) [27,47]. The mean *d*-excess value depended on the location and altitude of the station, but in all cases, higher *d*-excess monthly values were observed in autumn–winter precipitation, usually from September to December [47]. The lowest mean monthly *d*-excess values were observed in summer months. At the South Adriatic stations (Komiža, Dubrovnik), the distribution of the precipitation amount had higher seasonality (a low amount in summer and high in winter) than it did in Zagreb, and the *d*-excess values in summer months (from May to August) were lower than 8%, which is indicative of secondary evaporation of raindrops falling in a warm and dry atmosphere [47,50,51].

The very striking behavior of *d*-excess in March of the 2012–2018 period (Figure 11) can be explained as follows. During this period, the precipitation amount in March (57%) in four out of seven months was below 23 mm, resulting in *d*-excess values between −4.8% (2012, *P* = 4 mm) and 4.1% (2014, *P* = 22 mm), which was fully in accordance with observations from the station in Ljubljana [71] *d* < 5%) and corresponded with months with low precipitation. These *d*-excess values were low, but were still within 3σ of the mean value, and therefore they were not excluded from the analyses. When these lower *d*-excess values were excluded from the mean *d*-values for March from 2012 to 2018, the mean *d*-excess for 2012–2018 became 7.4%, which was not different from other periods. In contrast, in the 1980–2006 period, the total number of months with precipitation below 23 mm was only two (i.e., the occurrence of low precipitation in March was only 8%), and both *d*-excess values were excluded from the analyses. This observation also corroborated observed changes in the seasonal distribution of precipitation at continental stations in Croatia [59].

#### *4.4. Trends in Tritium Activity Concentration*

The seasonal variations in the tritium activity concentration in precipitation at Zagreb for the period from 1976 to 1993 were superposed onto a generally decreasing trend (Figure 5) that could be approximated by an exponential decay curve with a half-life of about 6 years, in accordance with the estimated residence time of tritiated water vapor in the lower stratosphere (of the order of a few years) [1]. The same pattern for the tritium-in-precipitation regime (maxima from June to July and minima in winter) was observed in other continental stations (Figure 1): Ljubljana [47,70,71], Plitvice Lakes [65,67], and Gacka [68]. The ratio of the maximum to minimum value for Zagreb precipitation up to 1996 ranged from 2.2 to 5.7 without a significant trend [62], similarly to other Northern Hemisphere stations (between 2.5 and 6 [1]). The ratio was similar in the 1996–2018 period, between 2.3 and 5.3 in most cases, and was higher only in years with tritium activity concentrations close to the GPC detection limit.

The decrease in mean annual tritium activity concentration values continued after 1996 (Figure 5 insert, Figure 12), but to a much lesser extent of 0.08 TU per year, resulting in mean values of 8.8 ± 1.4 TU in 1996–2006 and 7.6 ± 0.8 TU in 2012–2018. The mean values for the station Ljubljana were 9.1 TU and 8.3 TU during 1998–2010 and 2007–2010, respectively [71], showing the same trend.

**Figure 12.** Tritium activity concentration in precipitation at Zagreb during the period 1996–2018. The boxplot shows the median value and the percentile range 25%–75% in a shadowed box, — shows the percentile range 5%–95%, and the average value is shown by the symbol -.

Bomb-produced tritium in precipitation until about 1995 prevented studies on whether the natural production of tritium was influenced by variations in solar activities. The modulation of cosmogenic tritium production by an 11-year solar cycle has been recently shown in precipitation at several stations worldwide [93]. Local maxima in the tritium activity concentration in precipitation were observed simultaneously with maxima in neutron flux (minima in sunspot numbers). Our data (Figure 12) also showed local maxima in mean annual values and larger variability in 1996, 2007, and 2018 in accordance with the observations presented in Reference [93].

#### *4.5. Local Meteoric Water Line*

The slopes (*a*) and intercepts (*b*) of LMWLs were obtained using different regression methods (Table 3). The whole 1980–2018 period was taken, as were individual subperiods, and the present values were compared to the available data on LMWLs for Zagreb precipitation from different earlier periods (obtained using different regression methods). The slopes and intercepts determined using different regression methods increased from the OLSR to the RMA and MA (in the whole period, as well as in the subperiods (Table 3)), as was also observed for most continental stations [48]. The same was also valid for precipitation-weighted regressions. No difference was observed between the corresponding nonweighted and weighted types of regression, as was expected from the rather homogeneous distribution of the monthly precipitation amount. If the *rmSSEav* value was taken into account, all values were close to 1 (closer to 1 for RMA and PWRMA over the whole long-term period). In subperiod 1980–1985, PWMA represented the LMWL equally as well as PWRMA and RMA did; in subperiods 1986–1995 and 1996–2006, the best fits were obtained by RMA and PWMA; and in the most recent subperiods, the best fits were obtained by RMA and PWLSR (Table 3). It may be concluded that the local meteoric water line for Zagreb was best described by the nonweighted RMA regression method (Figure 13):

$$
\delta^2 \mathcal{H}\_{\text{RMA}, \text{all}} = (7.74 \pm 0.06) \,\delta^{18} \mathcal{O} + (5.6 \pm 0.6), n = 389. \tag{5}
$$

**Figure 13.** Local meteoric water line (LMWL) for Zagreb precipitation. Data from two periods are shown by different symbols. The fitted line is expressed by Equation (5).

A closer look at previously published LMWLs points to the difference in data up to 1996, with slope and intercept values close to 8 [27,62], while in 2001–2003, the slope was lower (7.3 + 0.2) and the intercept much lower (2.8 + 1.8) [47]. A similar conclusion was obtained from the present calculations (Table 3): all regression methods resulted in a slope close to 8 and an intercept in the range 7.3 to 9.3 in the subperiods 1980–1985 and 1986–1995. However, the slope values ranged from 7.4 to 7.8 and the intercept values from 2.6 to 6.3 in the subperiods 1996–2006 and 2012–2018. This difference can be explained by increases in temperature, especially in the summer months, and higher variability in the precipitation amount, which both led to more precipitation with lower *d*-excess, i.e., below the LMWL, as can be seen from Figure 13.

**Table 3.** Slopes (*a*) and intercepts (*b*) of the local meteoric water line (LMWL) for Zagreb in different periods and obtained using different regression methods. 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 datapoints included; *r* and *r*2: regression coefficients; *rmSSEav*: average of the root mean square sum of squared errors [48].


#### *4.6. Temperature Dependence of* δ*18O*

The relation between all values of the mean monthly air temperature *T* and the monthly δ18O values (Figure 14) can be described as

$$
\delta^{18}\text{O} = \left(0.331 \pm 0.013\right)T - \left(12.8 \pm 0.2\right),
\newline \text{ $n$ } = 394, \newline \text{ $r$ } = 0.795.\newline \tag{6}
$$

**Figure 14.** Temperature dependence of δ18O in monthly precipitation at the Zagreb station during different periods. Regression lines are represented by Equations (6) and (8). The mean values in the subperiods (stars) are taken from Table 2.

When data for the periods 1980–2006 and 2012–2018 were separated, the two relations were

$$
\delta \, ^{18}O\_{1980 \cdot 2006} = \left( 0.336 \pm 0.013 \right) T - \left( 12.8 \pm 0.2 \right), n = 320, r = 0.814,\tag{7}
$$

$$
\delta \, ^{18}O\_{2012\text{-}2018} = \text{(0.312} \pm 0.036\text{)} \, T - \text{(12.7} \pm 0.5\text{)}, \\
n = 74, \, r = 0.718. \tag{8}
$$

Although Equations (7) and (8) indicate a slight change in the temperature coefficient, the difference is comparable to the statistical uncertainties in the *a*-values, and therefore it can be concluded that there was no change in the relation δ18O versus *T* in the recent period compared to the older one. The two lines (Equations (7) and (8)) are practically indistinguishable (Figure 14). The new relation for the complete set of long-term data for Zagreb (Equation (6)) was not different (i.e., it is within uncertainties) from the same relation for the data for the 1980–1996 period (δ18O = (0.325 <sup>±</sup> 0.016) *T* − (12.6 ± 0.2), *r* = 0.83, *n* =183) [62]. The slopes in Equations (7) and (8) were also comparable to the slope of 0.31% per ◦C, which was determined from long-term data from midlatitude stations in the Northern Hemisphere [4].

#### **5. Conclusions**

A 43-year-long record of data on the isotope composition of precipitation (δ18O, δ2H, *d*-excess, and the tritium activity concentration *A*), together with meteorological data (air temperature, precipitation amount) at a continental station in Zagreb, Croatia, was studied, divided into four almost equally long subperiods. However, due to some missing data on stable isotope composition, the first (1980–1985) and the last (2012–2018) subperiods were shorter than the two subperiods in between (1986–1996, 1996–2006).

A constant increase in mean annual temperature was observed at a rate of 0.07 ◦C per year. An increase in monthly mean temperature was observed in all months in 2012–2018 when compared to earlier subperiods, with larger differences in the spring–summer than in autumn. The most striking feature of the annual precipitation amount was larger variations in the last two subperiods compared to earlier. A shift of the month with the highest precipitation amount was observed, from June and August to September.

Annual mean δ18O and δ2H values in the whole long-term period showed an increase of 0.017% per year and 0.14% per year, respectively. When monthly mean values in different subperiods were compared, larger differences were observed in the first half of the year than in the second one. Both of these changes in the stable isotope composition resembled observed changes in air temperature.

Although mean annual *d*-excess remained constant over the whole long-term period, there was a tendency for a decrease in the *d*-excess value in the first half of the year and an increase in the second half due to the influence of air masses originating from the eastern Mediterranean. Together with a shift in the maximal monthly mean value from October to November and a significant increase in *d*-excess in November, these observations point to changes in the precipitation regime and circulation pattern of air masses during the most recent period.

Different regression methods for the calculation of the local meteoric water line for Zagreb gave very similar values for the slope and intercept, with a slight preference for the RMA method and with no difference between nonweighted and precipitation-weighted values. In addition, no significant difference in both LMWLs and the temperature dependence of δ18O values was observed between the most recent period (2012–2018) and the earlier period (1980–2006). The observed temperature gradient of 0.33% per ◦C was comparable to that of other similar stations.

The tritium activity concentration in precipitation in Zagreb between 1976 and 1994 exhibited pronounced seasonal variations superposed on a generally decreasing trend with a half-life of about 6 years, which is typical for continental stations of the Northern Hemisphere. Since 1996, the mean annual *A* values have been almost constant, and the mean *A* value during the 2012–2018 period was 7.6 ± 0.8 TU. The tritium activity concentration in precipitation with no bomb-peak influence is worth further monitoring because of possible local contamination with technogenic tritium and also to enable studies of the solar cycle influence on the production of cosmogenic tritium.

The present analysis of long-term data on the isotope composition of precipitation may be useful for future comparisons to some other long-term records in nearby countries to obtain better knowledge on spatial and temporal variations across the wider region. It can also be a good basis for a comparison to some short-term records at other stations in Croatia. Last but not least, this analysis shows that climate changes are reflected in isotope compositions of precipitation, which means that further monitoring at stations with long-term records could be useful in studying the impact of climate changes on the environment, especially on water resources.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/12/1/226/s1, Figure S1: Monthly δ2H in precipitation at Zagreb, 1980–2018 period; Table S1: Stations from Croatia with available data on isotopes in precipitation; Table S2: Long-term monthly data (precipitation amount *P*, temperature *T*, δ18O, δ2H, *d*-excess, and tritium activity concentration *A*) for the station at Zagreb, Croatia; Table S3: Mean annual values δ18O, δ2H, *d*-excess, and *A*. Comparison of arithmetic and precipitation-weighted mean values.

**Author Contributions:** Conceptualization, I.K.B. and P.V.; data validation and interpretation, writing the manuscript—original draft preparation, final manuscript, visualization—most graphs, and communicating with the journal, I.K.B.; sample collection, tritium activity sample preparation, and discussion, J.B. and A.S.; data evaluation, database management, and LMWL Freeware software application, D.B.; stable isotope measurements, data reduction, and interpretation, P.V.; meteorological data analysis and visualization, I.L.M. All coauthors

*Water* **2020**, *12*, 226

participated in discussions and reviewed the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding. In the past, the isotope composition of precipitation was part of the following funders and projects: the Ministry of Science and Education (MSE) of the Republic of Croatia, project nos. 00980207 (natural radioisotopes and processes in gases, 1996–2002), 0098014 (natural isotopes of low-level activities and the development of instrumentation, 2002–2006), 098-0982709-2741 (natural radioisotopes in the investigation of Karst ecosystems and dating, 2007–2011); the IAEA, project nos. 11265 (Isotopic Composition of Precipitation in the Mediterranean Basin in Relation to Air Circulation Patterns and Climate Tritium and Stable Isotope Distribution in the Atmosphere in the Coastal Region of Croatia, 2000–2003), CRO/8/006 (the application of isotope techniques in the investigation of water resources and water protection in the Karst area of Croatia, 2005–2006 ), RER/8/012 (isotope methods for the management of drinking water resources in water-scarce areas, 2007–2008), CRO/8/007 (Using Isotope Tracers as a Tool for a Groundwater Vulnerability Assessment in the County of Split, Dalmatia, 2007–2008), RER/8/016 (Using Environmental Isotopes for an Evaluation of Streamwater/Groundwater Interactions in Selected Aquifers in the Danube Basin, 2009–2011), CRO/7/001 (Isotope Investigation of the Groundwater–Surface Water Interactions at the Well-Field Kosnica in the Area of the City of Zagreb, 2016–2017); EU project FP5 ICA2-CT-2002-10009 ANTHROPOL.PROT (a study of anthropogenic influence after the war and establishing protection measures in the National Park Plitvice and the Biha´c Region at the border area between Croatia and Bosnia and Herzegovina, 2003–2005); and a project financed by the National Park Plitvice Lakes ("an investigation of the influence of forest ecosystems of the National Park Plitvice Lakes on the quality of water and lakes", 2003–2006). Cooperation between Croatian and Slovenian partners was funded in recent decades by MSE and the Slovenian Research Agency (SRA) as part of bilateral cooperation (BI-HR/01-03-011, BI-HR/04-05-13, and BI-HR/09-10-032) and the SRA research program P1-0143.

**Acknowledgments:** The authors are thankful to former Ruđer Boškovi´c Institute laboratory staff members Bogomil Obeli´c, Nada Horvatinˇci´c, Dušan Srdoˇc, Elvira Hernaus, and Božica Mustaˇc and the present technician Anita Rajtari´c, who all took part in sample and data collection. We are also thankful to colleagues from the Jožef Stefan Institute in Ljubljana (Jože Pezdiˇc, Sonja Lojen, Nives Ogrinc, Silva Perko, Zdenka Trkov, and Stojan Žigon) for stable isotope measurements from the period 1980–2003 and to Jelena Parlov and Zoran Kovaˇc from the Laboratory for Spectroscopy of the Faculty of Mining, Geology, and Petroleum Engineering, University of Zagreb, for stable isotope measurements from the period 2012–2018. Meteorological data were obtained on request (free of charge) from the Croatian Meteorological and Hydrological Service. We thank Catherine E. Hughes and Jagoda Crawford for help in getting LMWL Freeware software and for discussions of the formulae used.

**Conflicts of Interest:** The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection; analyses; or interpretation of data; in the writing of the manuscript; or in the decision to publish the results

#### **List of Abbreviations**


#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
