Impact of Hydrological Conditions on the Isotopic Composition of the Sava River in the Area of the Zagreb Aquifer

: The Zagreb aquifer is the main source of potable water for the inhabitants of the City of Zagreb and Zagreb County. It presents a strategic water reserve protected by the Republic of Croatia. All previous studies related to the deﬁnition of the groundwater–surface interaction in the study area have been made based on the isotopic composition of the Sava River from the location of the Domovinski Most bridge, which is located downstream of most pumping well ﬁelds. In 2019, a new monitoring station was established at the Podsusedski Most bridge, at the entrance of the Sava River into the Zagreb aquifer, approximately 23 km upstream of the Domovinski Most bridge. Within this research, water isotope data ( δ 2 H, δ 18 O, deuterium excess) from both Sava River and groundwater sites were used along with hydrologic data to examine the extent to which hydrologic conditions affect the isotope signature and whether the interaction between groundwater and the Sava River causes a change in the isotopic composition of the Sava River. In addition, δ 18 O amplitudes were estimated for different time periods, as well as the mean residence time for the hydrological year 2019/2020. For that purpose, different statistical methods were applied to the new monthly data for six years for the Domovinski Most bridge and two years for the Podsusedski Most bridge. The δ 18 O amplitudes vary from 0.22 to 1.86 depending on the time interval and hydrological conditions, while the mean residence time for the hydrological year 2019/2020 was estimated to be about 2.5 months.


Introduction
Stable isotopes of hydrogen and oxygen from water are often used as a tracer for different hydrological processes. In the last decades they have been increasingly used due to advances in measurement techniques. They have become a very important tool in numerous hydrological research aims, including: the estimation of mean residence time (MRT) and groundwater dynamics [1][2][3][4][5], groundwater-surface water interaction [6,7], evaluation of surface water dynamics [8][9][10][11], evaluation of precipitation patterns [12][13][14][15][16][17], exploration of unsaturated zone water dynamics [18][19][20], and evaluation of other processes, primarily evaporation [2,[21][22][23]. In general, unconfined alluvial aquifers, which have connections with rivers, are dependent on the surface water fluctuations. It has been shown that many alluvial aquifers have a problem with groundwater depletion [24], while recently, that problem has been much more pronounced in the moderate climate areas, generally due to the impact of climate change [25,26]. Furthermore, it is well known that the relationship between groundwater and surface water can be one-way (gaining or losing stream), but it can also be very dynamic, depending mostly on the hydrological conditions, which can have big influence on changes in groundwater flow velocities and groundwater flow directions. These rapid changes make sometimes relevant hydrological relationships very hard to define. This highlights the necessity of studying the isotopic composition of the main recharge source in much more detail, especially if hydrologic relationships need to be explored in a local scale.
The Zagreb aquifer is designated as a part of a strategic water reserve, protected by the Republic of Croatia, which represents the main source of potable water for the citizens of the City of Zagreb and part of Zagreb County. In recent years, in the wider area of the City of Zagreb, numerous research projects have been performed, primarily isotopic, mostly focused on the following: the evaluation of nitrates' origin in the Zagreb aquifer [27], the influence of the Sava River temperature on the groundwater temperature [28], groundwater-surface water interaction between the Sava River and the Zagreb aquifer [6], the estimation of groundwater velocities [29], the evolution of Zagreb's Local Meteoric Water Line (LMWL), and changes in meteorological variables [16], as well as the examination of hillslope soil hydrology [30]. It has been shown that soil type can have big influence on precipitation infiltration [6], but also that a constant increase in mean annual air temperature is evident. Furthermore, in the last decades, annual precipitation amounts have varied much more than before [16]. There has also been research related to the evaluation of the Sava River water isotopic composition. In 2010, the Sava River isotopic signature was measured in Zagreb, near the Petruševec well field (Domovinski Most bridge), on a weekly time interval [31], as well as in 2016 at the same location, but on a monthly interval [6,32]. All this research showed that the Sava River was isotopically more like the precipitation which falls in Ljubljana, rather than in Zagreb, and that it has a very similar isotopic composition to the groundwater of the Zagreb aquifer. In the Sava River basin, recent research was related to the estimation of the spatial distribution of δ 18 O in precipitation, as well as to the distribution of oxygen and hydrogen stable isotopes and estimation of MRTs in Sava in Slovenia, upstream from the Zagreb aquifer, and in Serbia, in the area where the Sava River enters the Danube River [9]. Although it has been found that the main boundary condition related to the recharge of the Zagreb aquifer is the Sava River, which generates recharge from 67.5% up to 83.74% depending on the area which is evaluated [6], no research was focused on the evaluation of the Sava River isotopic pattern at the entrance of the Zagreb aquifer, i.e., the inflow area, and its isotopic difference with respect to the exit of the aquifer, i.e., the outflow area. The latest hydrogeological research confirmed that the hydraulic connection between the Sava River and groundwater is complicated and not uniform [33].
It has been proven that the Sava River is the main source of water for the Zagreb aquifer and influences the isotopic composition of the groundwater, but at the same time it can be assumed that the groundwater also influences the isotopic composition of the river, depending on the hydrological conditions. Despite numerous research efforts in the study area, the isotopic signature of the Sava River and its dependence on hydrologic conditions has not yet been studied in detail.
This research is focused on the evaluation of the isotopic composition of the Sava River at the entrance of the Zagreb aquifer, for which a new monitoring point has been established, as well as to its downstream part which can be approximated as the exit from the Zagreb aquifer. Furthermore, the main goals of this research are related to the calculation of the amplitudes of δ 18 O in the Sava River in different time intervals and different hydrological conditions, the comparison between groundwater and the Sava River isotopic signature, as well as to the estimation of the mean residence time between the two observed points in the Sava River.

Research Area
The Zagreb aquifer is the main source of drinking water for a quarter of the population of the Republic of Croatia. Due to its importance, it is protected by the Croatian state and designated as a strategic water reserve. It is located in the northwestern part of the Republic of Croatia and has six main well fields ( Figure 1). It consists of Quaternary sediments deposited in Middle and Upper Pleistocene and Holocene. The older sediments are mainly lacustrine and marshy deposits, while the Holocene deposits are alluvial. In general, the Pleistocene deposits are siliciclastic, while the Holocene are alluvial, which is due to the transportation of material from the Alps at the beginning of the Holocene, when the Sava River began to flow [34]. In the Holocene, gravels and sands dominate, while in the Pleistocene sediments, a frequent alternation of sands, gravels, silts and clays is observed [35].
From the hydrogeological point of view, the Zagreb aquifer is divided into the unsaturated part of the Zagreb aquifer, the shallow Holocene aquifer and the deeper Pleistocene aquifer.
The thickness of the unsaturated zone varies depending on hydrological conditions. Previous studies have shown that the thickness of the unsaturated zone varies between 2 and 13 m [6,36]. In addition, various soil types are found in the Zagreb aquifer area, with two main types dominating: Fluvisols and Eutric Cambisols on Holocene deposits [37,38]. Unfortunately, the unsaturated zone in a large part of the Zagreb aquifer is decomposed by human influence, especially on the left bank of the Sava River, where urban areas predominate.
The shallow Holocene aquifer is in direct contact with the Sava River, where the groundwater flow direction ranges from W/NW to E/SE. Previous studies have shown that the Sava River is the main recharge source for the Zagreb aquifer and that the influence of the Sava River is more pronounced in its vicinity [6,39]. It has also been shown that, in some areas, the Sava River drains the Zagreb aquifer at low and medium water It consists of Quaternary sediments deposited in Middle and Upper Pleistocene and Holocene. The older sediments are mainly lacustrine and marshy deposits, while the Holocene deposits are alluvial. In general, the Pleistocene deposits are siliciclastic, while the Holocene are alluvial, which is due to the transportation of material from the Alps at the beginning of the Holocene, when the Sava River began to flow [34]. In the Holocene, gravels and sands dominate, while in the Pleistocene sediments, a frequent alternation of sands, gravels, silts and clays is observed [35].
From the hydrogeological point of view, the Zagreb aquifer is divided into the unsaturated part of the Zagreb aquifer, the shallow Holocene aquifer and the deeper Pleistocene aquifer.
The thickness of the unsaturated zone varies depending on hydrological conditions. Previous studies have shown that the thickness of the unsaturated zone varies between 2 and 13 m [6,36]. In addition, various soil types are found in the Zagreb aquifer area, with two main types dominating: Fluvisols and Eutric Cambisols on Holocene deposits [37,38]. Unfortunately, the unsaturated zone in a large part of the Zagreb aquifer is decomposed by human influence, especially on the left bank of the Sava River, where urban areas predominate.
The shallow Holocene aquifer is in direct contact with the Sava River, where the groundwater flow direction ranges from W/NW to E/SE. Previous studies have shown that the Sava River is the main recharge source for the Zagreb aquifer and that the influence of the Sava River is more pronounced in its vicinity [6,39]. It has also been shown that, in some areas, the Sava River drains the Zagreb aquifer at low and medium water levels, while it releases water to the aquifer at high water levels [6]. The hydraulic boundaries of the Zagreb aquifer were determined based on the study of equipotential maps under different hydrological conditions [40]. A no-flow boundary was established in the north of the Zagreb aquifer, an inflow boundary in the south and west, and an outflow boundary in the east. Hydraulic conductivities range up to 3000 m/day, while the thickness of the shallow aquifer layer varies between 5 and 40 m [41,42]. Moreover, aerobic conditions prevail in most parts of the aquifer [43].
The deeper, Pleistocene aquifer layer has a thickness of up to 60 m in the eastern part, which means that the maximum thickness of the entire system is about 100 m [41]. Although the deeper aquifer layer is hydraulically connected to the shallower aquifer, geochemical stratification can be seen along the depth. Higher sodium concentrations were found in the deeper aquifer, probably due to lower groundwater velocities and longer residence times, and hydrogeochemical CaMgNa-HCO 3 facies formed compared to the hydrogeochemical CaMg-HCO 3 facies found in the shallow aquifer layer [44].
Various problems related to groundwater quality and quantity were observed. Previous studies identified several main groups of contaminants [42], while more recent studies focused on determining the origin and trend of nitrate [27,45] and the risk of nitrate contamination [46]. In terms of groundwater quantity, it was found that current groundwater levels in the Zagreb aquifer are generally about 3 to 6 m lower than historical groundwater levels observed in the 1960s [47]. There are numerous possible reasons for the lowering of groundwater levels, ranging from occasional flooding to the presence of embankments built along the Sava River that stopped a certain percentage of infiltration, to gravel extraction from the Sava River and excessive pumping of potable water. However, the most important reason is probably the extensive erosion of the riverbed caused by the regulation of the Sava upstream. In addition, climate change probably also contributes negatively to this problem. In recent decades, a steady increase in annual mean temperature and larger fluctuations in precipitation have been observed [16]. All this requires a detailed characterization of all variables related to the recharge component of the Zagreb aquifer if sustainable groundwater management is to be achieved.

Materials and Methods
Within this research, different hydrogeochemical data, as well as different statistical methods, have been used to evaluate the behavior of the Sava River in the inflow and outflow area of the Zagreb aquifer. Although the focus was on the evaluation of the water isotope signature (δ 2 H, δ 18 O, deuterium excess), in-situ chemical parameters and river water levels have been also observed.
For that purpose, previously published data as well as new data was used. Part of the groundwater and Sava River (Domovinski Most bridge-DM) isotope data was presented in Kovač et al. [27] and Parlov et al. [6], mainly related to the period from November 2015 till January 2017. All other data presents new, unpublished data (δ 2 H, δ 18 O, water temperature, pH, dissolved oxygen, and electrical conductivity). For the Domovinski Most bridge, located downstream (Figure 1), data used within this research is from January 2015 till March 2021 (in total 65 sampling campaigns). The data series end in March because in April 2019 it was decided that a new monitoring point of the Sava River would be established, in the inflow area of the Zagreb aquifer (Podsusedski Most bridge-PM). For the new monitoring point, the first two years of data is presented and used for the interpretation (in total 24 sampling campaigns). It must be emphasized that between these two monitoring points, the weir of the Zagreb cogeneration plant TE-TO is situated. Previous research showed that TE-TO Zagreb has strong influence on the dynamic interaction between groundwater and surface water [47], and consequently on groundwater flow velocities and directions.
All data has been evaluated in six main steps ( Figure 2). In the first step, basic statistical parameters (average, median, minimum, maximum, and standard deviation) related to the in-situ measurements have been observed (measured with WTW multi parameter 3630 IDS). In the second step, Sava water level frequency and duration curves were examined for the hydrologic station Podsused-Žičara which is located very close to the Podsusedski Most bridge. Within this step, frequency and duration curves from four hydrologic years (2016/2017, 2017/2018, 2018/2019, and 2019/2020) and for the longterm period (10/2000-9/2020) have been evaluated. Data for the Sava River water levels was provided by Croatian Meteorological and Hydrological Service. In the third step, water isotope data (δ 2 H, δ 18 O, deuterium excess) for both monitoring points at the Sava River, as well as for groundwater from previous studies, are presented through basic statistic parameters and bivariate plot. In the fourth step, boxplots have been created for water isotope data, after which extremes and outliers were firstly excluded from further analysis. It was assumed that the ejection of outliers and extreme values will generate more reliable results, i.e., better fit of seasonal sine wave curves. However, in the end, we decided to calculate amplitudes with and without outliers and extreme values. In the fifth step, data related to the Sava River and groundwater was interpretated quantitatively using periodic regression analysis to fit seasonal sine wave curves to annual δ 18 O variations as [9,48,49]:

Results and Discussion
In Table 1, basic statistical values of in-situ parameters of the Sava River at the two monitoring points are presented. Although the calculated values are similar, some differences can be seen. The temperature of the Sava River is slightly lower at the PM bridge, as well as values of electrical conductivity. Additionally, the standard deviation of all observed parameters is higher at the DM bridge. Higher values of temperature and electrical conductivity are expected downstream due to fact that Sava River flows through the urban part of the City of Zagreb. Furthermore, the average and median values of dissolved oxygen are lower at the DM bridge, although not significantly, which nonetheless indicates slightly higher mixing with groundwater. Previous research has shown that in the downstream part (i.e., eastern part), the Zagreb aquifer deepens, and dissolved oxygen concentrations are lower [42,43].  Estimation of the mean residence time (MRT) between the Sava River at the entrance (Podsusedski Most bridge) and exit (Domovinski Most bridge) of the Zagreb aquifer was performed in the sixth step. For the calculation of the MRT, a commonly used exponential model was used [9,[48][49][50]: where A PM is the fitted amplitude in the area of the PM bridge, A DM is the fitted amplitude in the area of the DM bridge, while c is the radial frequency of the annual fluctuation as presented in Equation (1) Figure 1 was constructed using ArcMap 10.8.1., while georeferenced orthophoto background was obtained from the geoportal of the Croatian Geodetic Administration (map is presented using the official coordinate system of the Republic of Croatia-HTRS96/TM).

Results and Discussion
In Table 1, basic statistical values of in-situ parameters of the Sava River at the two monitoring points are presented. Although the calculated values are similar, some differences can be seen. The temperature of the Sava River is slightly lower at the PM bridge, as well as values of electrical conductivity. Additionally, the standard deviation of all observed parameters is higher at the DM bridge. Higher values of temperature and electrical conductivity are expected downstream due to fact that Sava River flows through the urban part of the City of Zagreb. Furthermore, the average and median values of dissolved oxygen are lower at the DM bridge, although not significantly, which nonetheless indicates slightly higher mixing with groundwater. Previous research has shown that in the downstream part (i.e., eastern part), the Zagreb aquifer deepens, and dissolved oxygen concentrations are lower [42,43].  Figure 3 shows the flow duration curve for the long-term data (2000 to 2020) and the last four subsequent hydrologic years. It is immediately apparent that the larger flows (greater than 5% of the year) exhibit quite a bit of variability compared to the long-term data curve. The mid flows (between 5% and 65% of the year) are lower in every year except 2017/2018, while lower flows (greater than 65% of the year) have been above the long-term average in all four years. Moreover, the hydrological year 2017/2018 represents an extremely wet year with a very high duration of almost all water levels of the Sava River. Figure 4 shows the frequency of water levels of the Sava River for four hydrological years.
The year 2017/2018 is clearly different from the others, which can be also seen in the much more frequent occurrence of the Sava River water levels reaching between 118 m a.s.l. and 121 m a.s.l. more frequent occurrence of the Sava River water levels reaching between 118 m a.s.l. and 121 m a.s.l.   Table 2, basic statistical parameters related to the isotopic composition of the groundwater of the Zagreb aquifer and the Sava River are shown. The results suggest that the groundwater isotopic signature is slightly different with respect to the Sava River.   In Table 2, basic statistical parameters related to the isotopic composition of the groundwater of the Zagreb aquifer and the Sava River are shown. The results suggest that the groundwater isotopic signature is slightly different with respect to the Sava River. In Table 2, basic statistical parameters related to the isotopic composition of the groundwater of the Zagreb aquifer and the Sava River are shown. The results suggest that the groundwater isotopic signature is slightly different with respect to the Sava River. However, the average and median values from the DM bridge are closer to the isotopic signature of the groundwater, while d-excess values are almost identical. This corresponds to the visual inspection that can be drawn from Figure 5, in which it can be also seen that the Sava River isotopic composition from the DM bridge is more similar to groundwater from the Zagreb aquifer. Furthermore, the results show much more deviation in the isotopic signature at the DM bridge. Together with the results of in-situ measurements and previous research, where it was shown that the connection between Sava River and Zagreb aquifer system is not uniform [6,33,42,43], these results indicate that mixing between surface water and groundwater is more pronounced in the downstream part of the Zagreb aquifer. Greater aquifer depth in the downstream part, as well as the existence of TE-TO Zagreb, generates slower velocities and enables more time for infiltration of the Sava River into the Zagreb aquifer, which is consistent with the results of the previous research [6]. Within this research, preliminary results related to the evaluation of the water isotope composition suggested the existence of outliers and extreme values. In order to quantify these and remove them from further analysis, boxplots were constructed. Groundwater has the greatest data set, and in Figure 6 it can be seen that it also has the most outliers and extreme values. In the next steps, only values of δ 18 O were used. The boxplots shown in Figure 6 suggest that only 3 values of δ 18 O for the DM bridge should be removed, 13 for groundwater, and none for the PM bridge. It is very interesting to see that all outliers and extreme values for the DM bridge are related to the time period of spring/summer 2018, i.e., part of the hydrological year 2017/2018, which was an extremely wet hydrological year, but also a very warm to extremely warm year, according to the average monthly temperature values as presented by the Croatian Meteorological and Hydrological Service (meteo.hr, accessed on 8 June 2022). The exclusion of these data would make it impossible to perform periodic regression analysis to fit seasonal sine wave curves to annual δ 18 O variations for the hydrologic year 2017/2018. In the end, it was decided to use outliers and extreme values which were observed in hydrological year 2017/2018, which were found to be very useful for the interpretation and showed that the extreme isotopic signature corresponds to the occurrence of an extremely wet and hot hydrological year. Within this research, preliminary results related to the evaluation of the water isotope composition suggested the existence of outliers and extreme values. In order to quantify these and remove them from further analysis, boxplots were constructed. Groundwater has the greatest data set, and in Figure 6 it can be seen that it also has the most outliers and extreme values. In the next steps, only values of δ 18 O were used. The boxplots shown in Figure 6 suggest that only 3 values of δ 18 O for the DM bridge should be removed, 13 for groundwater, and none for the PM bridge. It is very interesting to see that all outliers and extreme values for the DM bridge are related to the time period of spring/summer 2018, i.e., part of the hydrological year 2017/2018, which was an extremely wet hydrological year, but also a very warm to extremely warm year, according to the average monthly temperature values as presented by the Croatian Meteorological and Hydrological Service (meteo.hr, accessed on 8 June 2022). The exclusion of these data would make it impossible to perform periodic regression analysis to fit seasonal sine wave curves to annual δ 18 O variations for the hydrologic year 2017/2018. In the end, it was decided to use outliers and extreme values which were observed in hydrological year 2017/2018, which were found to be very useful for the interpretation and showed that the extreme isotopic signature corresponds to the occurrence of an extremely wet and hot hydrological year. All models are robust and statistically significant, with R 2 values generally higher than 0.5. The deviation is observed only for the extreme hydrological year 2017/2018 which has several times higher amplitude than all other hydrological years (1.86 with respect to 0.4, 0.22, and 0.24). However, it must be emphasized that this model is not as robust as the others, with p value close to 0.05 and R 2 of 0.4. Slightly higher amplitudes in hydrological year 2016/2017 are probably the consequence of the more frequent occurrence of the Sava River water levels rising above 120 m a.s.l. Periodic regression analysis of groundwater amplitudes of δ 18 O (not shown) resulted in very small amplitudes (0.02 without outliers and extremes and 0.07 with outliers and extremes) with statistically insignificant models which have very small R 2 values (0.05 and 0.21, respectively). Despite the fact that more years of data are needed for the evaluation of isotopic signature of precipitation, the results suggest that the most appropriate method for the evaluation of the isotopic signature of the Sava River is based on the hydrological year. It was shown that δ 18 O amplitudes can vary greatly due to different hydrological conditions. Furthermore, the results suggest that MRT for hydrological year 2019/2020 was about 2.5 months (approximately 80 days). However, continuation of isotope monitoring is necessity if estimation of the amplitude and MRT values want to be calculated in different hydrological conditions and with greater reliability. Periodic regression analysis was used to fit seasonal sine wave curves to annual δ 18 O variations. In Figure 7a,b, amplitudes for the PM bridge are shown. Figure 7a presents an

Conclusions
This research presents new findings related to the isotopic composition of the Sava River at the entrance and exit of the Zagreb aquifer, which is a very important aquifer protected by the Croatian state as a strategic water reserve. Old and new isotopic data, in situ parameters, and detailed hydrological and statistical analyses were used to identify the impact of hydrological conditions on the isotopic signature of the Sava River in the area of the Zagreb aquifer. In detail, considering the observational data, the following results are indicated:

•
The isotopic signature together with the observed in situ parameters indicate that mixing with groundwater is more pronounced in the downstream part of the Zagreb aquifer. • Evaluation of δ 18 O amplitudes in the Sava River showed that they vary between 0.22 and 1.86, mainly depending on hydrological conditions. The extreme isotopic signature was the result of the extremely wet and hot hydrological year 2017/2018. • Although long-term data can generate reliable results in the evaluation of the precipitation isotopic signature, it has been shown that the isotopic signature in the Sava River must be evaluated at the hydrologic year level and that the use of outliers and extreme values should be studied in detail.

•
The MRT for the 2019/2020 hydrologic year was estimated to be about 2.5 months.
Finally, it can be concluded that not only does the Sava River influence the isotopic signature of groundwater in the study area, but also that groundwater influences the isotopic signature of the Sava River under certain hydrological conditions. Author Contributions: Conceptualization, Z.K. and J.B.; methodology, Z.K.; formal analysis, Z.K. and J.P.; investigation, Z.K., J.B., J.P. and A.S.; data curation, Z.K. and J.P.; writing-original draft preparation, Z.K.; writing-review and editing, J.B., J.P. and A.S.; visualization, Z.K. and J.P.; supervision, J.B. and J.P.; sampling and analysis, Z.K., J.B., J.P. and A.S. All authors have read and agreed to the published version of the manuscript.