Application of Hydrochemical Parameters as Tool for Sustainable Management of Water Supply Network

Abstract

Effective management of Water Supply Systems (WSSs) is crucial for ensuring safe drinking water. The WSS of Varaždin County is a complex network involving three groundwater sources: Bartolovec and Vinokovščak wellfields (alluvial aquifers) and Bela karstic spring. To achieve a comprehensive characterization of WSSs, routine laboratory data was integrated with stable isotopes and geochemical modeling. Within this study, all measured parameters remain below the Maximum Contaminant Level (MCL), ensuring water safety for human consumption. The Piper diagram identified variations in water sources based on their chemical composition, providing a simplified overview of mixing patterns within WSSs. Among the modeling approaches, inverse modeling (IM) was found to be more reliable than forward modeling (FM) and mass balance modeling (MB). Despite the limited capacity of δ¹⁸O to provide accurate mixing results, it was revealed that the reservoir water was in equilibrium with the air (no evaporation effects), indicating well-sealed reservoirs. Mixing modeling showed that the western, southwestern, and northern parts of the WSS mixed all three sources, whereas the eastern and southeastern areas primarily relied on the deeper aquifer of the Bartolovec source, indicating potential vulnerability. Strict validation criteria ensured the reliability of results, demonstrating the effectiveness and applicability of geochemical modeling in water security management plans.

1. Introduction

Water for human consumption is a precious resource that affects the development of civilization. While often taken for granted in developed nations, access to safe drinking water remains a significant challenge in less developed or underdeveloped countries, particularly among marginalized communities. According to the World Health Organization [1], over 2 billion people do not have access to safe water for human consumption. Consequently, ensuring access to drinking water for the population becomes progressively more challenging and expensive due to growing demand, the distance of available water source(s) from urban areas, enhanced contamination risks, and climate change impacts [2]. To address these challenges, municipal water managers rely on numerous technologies to ensure the capacity to deliver water with adequate quality and quantity to their consumers [3]. Water suppliers monitor water to ensure its quality meets the standards according to the EU Directive 2020/2184 on the quality of water intended for human consumption [4]. Therefore, it is necessary to protect and manage not only water resources but also the WSS itself. Effective management of the WSS not only reduces water losses (ensuring quantity) but also protects water from contamination within the system.

Within the WSS, water passes through an intricate network of pipelines, valves, latches, and water reservoirs (storage tanks); before reaching end-users, various changes can affect the water quality and, consequently, customer satisfaction with the delivered water. Therefore, it is necessary to have a comprehensive understanding of the system to quickly address incidents such as pipeline ruptures, water contamination, or user dissatisfaction. This knowledge enables the identification of causes and the implementation of control measures to mitigate risks to acceptable levels. To evaluate water quantity within a water supply network, suppliers commonly employ deterministic or probabilistic analyses to determine the flow rates and pressures in different sections of the network [3,5,6]. Generally, such analysis is robust; it is computationally intensive due to the lack of unique solutions, and the accuracy of these analyses is difficult to verify due to the lack of field measurements within the networks [7].

The application of stable isotopes (δ¹⁸O, δ²H) and geochemical parameters have been widely used as natural tracers in delineating the hydrological cycle of water to better understand the hydrogeological systems [8,9,10,11,12]. In recent studies, the application of stable isotopes and geochemical parameters was used to determine the status of WSSs because they are methodologically very straightforward and allow data to be uniquely connected to the consumer [13,14,15,16]. These parameters were shown to be more reliable compared to simplified deterministic or probabilistic network analyses [3,7]. Furthermore, various types of modeling, including mixing models and both forward and inverse geochemical models utilizing isotopic and geochemical parameters, have been used to obtain information about leaks in the water supply network, sewage permeability, and even correlations with a socioeconomic indicator such as population size and income [13].

To the best of our knowledge, the approach of determining the water transport and mixing within the network has not been used in Croatia. Only in recent years has research been conducted to explore the effects of changing water sources on water quality in the drinking WSS [17,18]. However, the application of stable isotopes has been widely used to determine groundwater age, identify sources of contamination, analyze surface and groundwater interaction, etc. [19,20,21,22,23,24,25,26]. In addition, the validation of mixing models using selected chemical parameters measured in initial and end member waters helps to estimate the reliability of gained values by mixing modeling.

The aim of this paper is to employ a multi-proxy approach that integrates geochemical parameters (major ion chemical data and in-situ physiochemical data), stable isotopes (δ¹⁸O, δ²H), and different types of mixing models to analyze water flow paths within the network and to determine the origin of water in different sections of the network. In addition, this study compares different mixing models and provides recommendations on which model is best suited for complex systems such as in the Varaždin area. The major novelty of this work is demonstrating that simple tools, specifically geochemical parameters that are routinely measured in water supply laboratories, combined with modeling and stable isotopes, can be useful tools for the management of water supply networks, not just on a local but also on a global scale. Furthermore, this study introduces a novel method for validating mixing models, enhancing the reliability of their results.

2. Materials and Methods

2.1. Study Area

This study focuses on the Varaždin area in NW Croatia, where potable water is distributed by the public water supplier Varkom d.o.o. (Figure 1). The main source of potable water is groundwater extracted from the Varaždin aquifer at two wellfield sites, Bartolovec and Vinokovščak, as well as a spring in Belski dol situated at the foothill of Mt. Ivanščica. Water from these three sources is, after disinfection with chlorine, transported to consumers via a 1670 km long distribution network.

The Bartolovec wellfield site, situated approximately 7 km east of Varaždin town in the village of Bartolovec, serves as the primary water source for the majority of consumers in Varaždin County, meeting around 70% of the water needs for human consumption in the regional WSS. The annual average rate of pumped water is 264 L/s (data for 2022) from 9 wells (B1-B9), of which 5 capture water from the upper aquifer (depth up to 45 m) and 4 from the lower aquifer (depth up to 120 m). The Quaternary-age alluvial aquifer in this area is divided by a 5 m thick silty-clay aquitard between the shallow aquifer (SA) and the deeper aquifer (DA) [27]. This aquitard protects the DA from pollution from the surface due to its poor permeability. Pumped groundwater from wells is disinfected and mixed in a ring-shaped pipeline system with five exits designated for distribution: Varaždin, Ludbreg, Novakovec, Doljan, and Tonimir (Figure 1). Because water mixes in a ring-shaped pipeline, it is hard to know the specific source and quantity of water distributed in each direction. This information is especially important in case of incidents when potential deterioration of water quality in the water supply network could occur.

The Vinokovščak wellfield site is located approximately 3 km NW of Varaždin, where groundwater is pumped from 4 wells (ZV1-ZV4), with an average annual rate of 82 L/s (data for 2022, Figure 1). Due to the thinner alluvial aquifer in this area, the wells draw water from a maximum depth of 45 m. Screens in wells are installed along the depth of wells, facilitating the mixing of water from SA and DA during pumping. After disinfection, the pumped groundwater from wells is distributed to one pipeline towards the north part of Varaždin town and the NW part of Varaždin County.

Spring water from Bela is captured at an average rate of 30 L/s (Figure 1). It is located in the village of Bela at the foot of the Ivanščica Mountain, 20 km SW from the town of Varaždin. This water source gravitationally flows through the pipeline to the pumping station Filipići, where it is disinfected and distributed into the water supply network. The Bela spring represents the smallest groundwater source within the WSS, covering 8% of the region’s water needs. Unlike the previously described groundwater sources, the catchment area of the spring Bela consists of carbonate-fissured rocks, mainly dolomite of the Triassic age [28].

2.2. Water Sampling

Raw water samples (prior disinfection) were collected on a monthly basis between December 2021 and December 2022 from groundwater sources at Bartolovec, Vinokovščak, and Bela for chemical and isotopic analyses (Figure 1, Table 1). At the Vinokovščak wellfield, one well (ZV-4) was sampled, while at the Bartolovec wellfield, samples were taken from three wells: two from SA (B-1 and Zp-7) and one from DA (Zd-4). Additionally, at the Bartolovec wellfield, samples were collected from all five (Figure 1, Table 1). Sampling from the distribution network was carried out twice, in July and November 2022. During the July 2022 campaign, reservoirs and public taps were sampled, while in the November 2022 campaign, reservoirs and hydrants were sampled (Figure 1, Table 1). Since the second campaign occurred during winter, public taps were closed to prevent pipe bursting due to low air temperatures.

Samples were poured into a 50 mL plastic bottle with a tight-fitting cap for stable isotope analyses and in 350 mL for chemical analyses. All samples were measured in the laboratory immediately upon returning from the field.

2.3. Water Analyses

Field measurements of electrical conductivity (EC), water temperature (T), pH, and dissolved oxygen content (DO) were conducted using a calibrated portable WTW instrument. Measured DO values for wells were excluded from consideration because they were affected by pumping. Total alkalinity was determined in the field by titration with 1.6 N H₂SO₄ using phenolphthalein and bromocresol green-methyl red as indicators. Major cations and anions were analyzed using the ion chromatography technique on a Dionex ICS-6000 DC at the Hydrochemical Laboratory of the Croatian Geological Survey. The analytical precision of cation and anion measurements, indicated by the ionic balance error (IBE), was computed on the basis of ions expressed in meq/l. The IBE value was observed to be within a limit of ±5% [29,30].

The δ¹⁸O and δ²H were determined using the Picarro L2130i (Santa, Clara, USA) in the Hydrochemical Laboratory of the Croatian Geological Survey. The instrument uses CRDS (Cavity Ring-Down Spectroscopy) technology [31]. All measurements were checked against USGS standards (USGS48 −2.224 ± 0.012‰ δ¹⁸O; −2.0 ± 0.4‰ δ²H; USGS46a −30.09 ± 0.05‰ δ¹⁸O; −235.6 ± 0.7‰ δ²H), which were checked periodically against the International Atomic Energy Agency (IAEA) standards: Vienna Standard Mean Ocean Water 2 (VSMOW2) and Standard Light Antarctic Precipitation 2 (SLAP2). Measurement precision was ±0.2‰ for δ¹⁸O and ±1‰ for δ²H. It is generally known that all isotopic results are expressed as per the international measurement standard, VSMOW2 [32,33].

2.4. Modeling

This study employed three approaches to calculate water mixing ratios: the mixing models. The first approach involves the application of mass balance calculations (MB), a commonly used method [34,35,36] where proportions of two-end members in the observed sample can be calculated using simple formulas based on conservative tracers Cl⁻ and δ¹⁸O.

Specifically, the calculations were performed using the following two equations:

(1)$f_{S1}={\displaystyle \frac{{\mathrm{C}\mathrm{l}}_{\mathrm{s}\mathrm{a}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{e}}^{-}{-\mathrm{C}\mathrm{l}}_{\mathrm{s}2}^{-}}{{\mathrm{C}\mathrm{l}}_{\mathrm{s}1}^{-}{-\mathrm{C}\mathrm{l}}_{\mathrm{s}2}^{-}}}$

(2)$f_{S1}={\displaystyle \frac{{\mathsf{\delta }}^{18}{\mathrm{O}}_{\mathrm{s}\mathrm{a}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{e}}{-\mathsf{\delta }}^{18}{\mathrm{O}}_{\mathrm{s}2}}{{\mathsf{\delta }}^{18}{\mathrm{O}}_{\mathrm{s}1}{-\mathsf{\delta }}^{18}{\mathrm{O}}_{\mathrm{s}2}}}$

where f_S1 represents the fraction (between 0 and 1) of source water 1 estimated in a sample of mixed origin. The Cl⁻_sample and δ¹⁸O_sample represent concentrations in waters from reservoirs, public taps, exits, and hydrants; Cl⁻_s1, δ¹⁸O_s1, Cl⁻_s2, and δ¹⁸O_s2 represent concentrations in different groundwater sources.

The other two models were conducted using geochemical modeling software, which uses geochemical parameters to infer the physical and chemical processes that control water chemistry along flow paths, specifically within the water supply network (pipeline, reservoirs), and to quantify reactions occurring within this network. The inverse geochemical modeling (IM) was performed using the PHREEQC software, version 3-A [37]. It is based on a geochemical mole-balance program that calculates the phase mole transfers (the moles of minerals and gases that enter or leave a solution) to account for the differences observed between initial and final water composition along the flow path, in our case, within the network. The third approach involved using the NETPATH-WIN software for forward geochemical modeling (FM) [38,39]. In all three models, EC, Cl⁻, NO₃⁻, and δ¹⁸O were used as constraints. Each of them has its advantages and disadvantages, which will be explained in more detail in the Results and Discussion chapter. The accuracy of the results was controlled according to the saturation index data for calcite and dolomite, which were compared with the actual values.

At the Bartolovec source, mixing models were performed in two steps because the mixing of SA and DA occurs at the source prior to water reaching the water supply network. The first step was to calculate the mixing proportions of groundwater from SA (represented by B-1 and Zp-7) and DA (represented by Zd-4) in the final waters of the wellfield, which are the exits (E1-E5). For the groundwater sources Vinokovšćak and Bela, this step was unnecessary due to the presence of only one exit into the network. In the second step, mixing models in the network (reservoirs, hydrants, and public taps) were calculated, where exits from all water sources (E1-E5, ZV-4, and Bela) represented initial waters, and samples from the network (H1-7, PT1-7, and R1-13) represented final waters. The modeling results are presented as percentage ratios of the initial waters in the final ones (ranging from 0 to 100%).

The validation of all three mixing models was performed by calculating sulfate concentrations in end members and comparing them with measured values. The sulfate concentration was chosen due to its conservative nature, inability to precipitate or to be added during disinfection, and it was not used as a tracer in calculations. Considering a 5% measurement error on IC, all values below 0.5 mg/L were taken as reliable results, while values exceeding this threshold were considered as poor results.

3. Results and Discussion

All measured parameters in raw water (groundwater sources) and water supply network, including mean (average), minimum, maximum, and standard deviation, are presented in Tables S1 and S2 in the Supplementary Material.

3.1. Geochemical and Isotopic Characteristics of Groundwater Sources

3.1.1. Bartolovec Wellfield

Measured temperature values in SA varied from 12.5–14.9 °C for well B-1 and 12.2–13.3 °C for well Zp-7 (Table S1). Although both wells capture the same aquifer, the observed difference in temperatures is attributable to the shallower depth of well B-1 (32 m) compared to well Zp-7 (45 m). Consequently, well B-1 is susceptible to water percolation from the surface during different seasons, while deeper Zp-7 remains unaffected by such phenomena. The average measured temperature in well Zd-4 during 12 months was 13.2 °C, indicating homogenization and longer transit time of water within the DA. Measured pH values showed differences between SA and DA. The pH values in SA varied from 7.2 to 7.7, showing higher oscillations throughout the year, while pH values in DA are slightly higher, ranging from 7.4 to 7.7 without significant oscillations during the year (Table S1).

Measured EC and nitrate concentrations in SA showed higher oscillations throughout the year compared to DA (Figure 2). Nitrate concentrations in well B-1 varied from 13.7 to 22.4 mg/L and in well Zp-7 from 19.0 to 26.0 mg/L, demonstrating significant oscillations during the year (Figure 2). The intensive agricultural production in the catchment area of these wells, coupled with a very thin surface-covering layer, contributes to leaching nitrates from the surface into the groundwater. However, this phenomenon is not observed in the well Zd-4, where nitrate concentrations range from 4.7 to 9.2 mg/L, additionally indicating a longer transit time of water in DA compared to SA. This observation confirms that the clay-silty layer that divides the aquifers has a protective function [25,40]. The measured EC values were higher in SA compared to DA. In SA (wells B-1 and Zp-7), EC values varied from 650 to 680 µS/cm, while in DA, EC varied from 469 to 480 µS/cm (Figure 2). This discrepancy is a consequence of higher concentration of dissolved solids in waters of SA relative to DA. In addition, the measured chloride concentrations in SA are significantly higher compared to DA. The average chloride concentration in well B-1 is 27.1 mg/L, and in well Zp-7, 24.7 mg/L, whereas in well Zd-4, the average is only 3.42 mg/L. High chloride and nitrate concentrations in the SA indicate washing out from the unsaturated zone, as previously observed in the study by Karlović et al. [41].

Geochemical differences between SA and DA were also visible in their Ca²⁺, Mg²⁺, SO₄²⁻, Na⁺, K⁺, HCO₃⁻, and F⁻ concentrations. In SA, these parameters are higher, with average values of 94.3 mg/L, 17.3 mg/L, 31.7 mg/L, 14.1 mg/L, 4.24 mg/L, 336.7 mg/L, and 0.39 mg/L, respectively, compared to DA, where the corresponding averages are 73.9 mg/L, 13.8 mg/L, 6.94 mg/L, 5.21 mg/L, 1.31 mg/L, 304.8 mg/L, and 0.30 mg/L (Table S1). According to the Piper diagram (Figure 3), both waters belong to the CaMg-HCO₃⁻ hydrochemical type, which is a consequence of the dissolution of carbonate and silicate minerals that form the aquifer [27,42,43]. However, even on this diagram, differences between waters are visible.

In addition to the aforementioned parameters, monthly measurements of stable oxygen and hydrogen isotopes in water also revealed differences between SA and DA, especially in the warmer part of the year (Figure 4). In SA, δ¹⁸O values range from −10.69 to −9.55 ‰ and δ²H from −74.5 to −67.8 ‰, while in DA, they range from δ¹⁸O −10.69 to −10.1 ‰ for δ¹⁸O and from −73.6 to −70.6 ‰ for δ²H (Figure 4, Table S1). Variations between wells B-1 and Zp-7, although they capture the same aquifer, mirror the patterns observed with temperature, chloride, and nitrate concentrations. This phenomenon underscores the heterogeneous nature of SA. All these differences provide valuable insights for developing mixing models.

3.1.2. Vinokovščak

The measured EC values of groundwater at the Vinokovščak wellfield ranged from 622 to 658 µS/cm (Figure 5). The temperature values ranged from 11.9 to 13.1 °C, while pH values ranged from 7.24 to 7.48 (Table S1). Chloride and nitrate concentrations ranged from 7.6 to 10.9 mg/L and 16.1 to 25.5 mg/L, respectively, without significant oscillations, although a decreasing trend was observed during the monitored period (Figure 5). The monitored year experienced very dry conditions, and this part of the study area is very well connected with the Drava River, especially during periods of low groundwater levels [41]. The Drava River has low concentrations of chloride and nitrate [41], which could contribute to the observed decreasing trend due to the higher river influence towards hydrological minimums.

The concentrations of major ions (Ca²⁺, Mg²⁺, SO₄²⁻, Na⁺, K⁺, and HCO₃⁻) also did not show significant oscillation, with measured average values of 98.5 mg/L, 19.3 mg/L, 29.2 mg/L, 4.14 mg/L, 1.52 mg/L, and 357 mg/L, respectively. According to the Piper diagram (Figure 3), monitored waters belong to the CaMg-HCO₃⁻ hydrochemical type, attributed to the dissolution of carbonate and silicate minerals that build the aquifer [27,42,43]. Interestingly, these waters have a similar Ca²⁺/Mg²⁺ ratio as waters from deep well Zd-4 of the Bartolovec wellfield, but chloride and sulfate concentrations are higher, reflecting different geochemical conditions within the aquifer (Figure 3).

On the contrary, measured stable oxygen and hydrogen values closely resemble those of the Bartolovec wellfield (Figure 6). It was observed that δ¹⁸O oscillated within the range of ±0.2 ‰, consistent with the measurement error. Similarly, δ²H varied around measurement error of ±1 ‰_. This initial observation was the first hint that δ¹⁸O would not be the most suitable parameter for mixing modeling. Furthermore, some values are plotted above the LMWL in the δ¹⁸O–δ²H diagram, especially wells DA and ZV-4, which is explained in detail in [11].

3.1.3. Bela

All measured hydrochemical parameters at this groundwater source are a consequence of the long transient time of the groundwater [44]. The measured water temperature ranged from 11.1 to 11.3 °C, while the pH values ranged from 7.4 to 7.6 (Table S1). The same behavior is observed in EC values, ranging from 495 to 501 µS/cm (Figure 7). Higher oscillations were observed in chloride, nitrate, and sulfate concentrations (Figure 7, Table S1), which is attributable to the influence of surface water. The spring is situated along the Varaždin–Podrute road, and it is divided by road by stream. During the high-water periods, surface runoff into the spring occurs.

Compared to the other two groundwater sources, the spring water contains high concentrations of Mg²⁺ (22.9–25.6 mg/L) due to the dolomite dissolution in the catchment area of the spring. Measured average values of other major ions Ca²⁺, Na⁺, K⁺, and HCO₃⁻ are 62.5 mg/L, 1.16 mg/L, 0,49 mg/L, and 314.4 mg/L, respectively. According to the Piper diagram, spring water belongs to the CaMg-HCO₃⁻ hydrochemical type (Figure 3), noticeably different compared to the other two sources.

The measured δ¹⁸O values range from −10.67 to −10.29 ‰, and the δ²H ranged from −72.20 to −70.80 ‰ (Table S1). These values are the most negative of all investigated groundwater sources and are shifted from the LMWL for Varaždin (Figure 6). This discrepancy can be attributed to several factors. Firstly, the catchment area of the spring is situated at a higher elevation (averaging about 700 m a.s.l.) with the highest peak exceeding 1000 m a.s.l., resulting in higher, more negative δ¹⁸O values in spring water. Secondly, the transit time in such aquifer types is long, over 50 years [44]. During this period, the winters were colder, with high amounts of snow [45] serving as the major source of recharge.

3.2. Geochemical and Isotopic Characteristics of Water from Water Supply Network—Reservoirs, Public Taps and Hydrants

3.2.1. Reservoirs

According to the measured data for reservoirs, pH values ranged from 7.3 to 8.2, temperatures from 11 to 21 °C, DO from 0.8 to 10.4 mg/L, EC from 486 to 683 µS/cm, and HCO₃⁻ from 295 to 360 mg/L. In reservoir water, high oscillation of basic cations and anions was observed due to the mixing of different groundwater sources. An important finding was that none of the measured concentrations exceeded those measured in the groundwater source, indicating no significant impact (contamination) within the network. Measured values ranged as follows: Ca²⁺ from 61.3 to 96.8 mg/L; Mg²⁺ from 13.7 to 26.8 mg/L; Na⁺ from 1.3 to 14.9 mg/L; K⁺ from 0.4 to 4.47 mg/L; Cl⁻ from 2.94 to 22.4 mg/L; SO₄²⁻ from 8.86 to 32.4 mg/L; and NO₃⁻ from 2.5 to 22.4 mg/L. The measured stable isotope values for δ¹⁸O varied from −10.67 to −9.81 ‰, and for δ²H from −71.9 to −69.3 ‰.

3.2.2. Public Taps and Hydrants

Similar to the reservoirs, the geochemical parameters of sampled waters from public taps and hydrants showed variability, but concentrations did not exceed the maximum values of groundwater sources. The measured pH values ranged from 7.28 to 7.69, EC from 483 to 676 µS/cm, temperature from 8.3 to 23.1 °C, DO from 0.8 to 8.9 mg/L, and HCO₃⁻ from 302 to 366 mg/L. The concentrations of major ions were as follows: Ca²⁺ from 71.8 to 98.2 mg/L; Mg²⁺ from 13.8 to 19.0 mg/L; Na⁺ from 4.0 to 15.0 mg/L, K⁺ from 1.4 to 4.5 mg/L; Cl⁻ from 4.5 to 22.8 mg/L; SO₄²⁻ from 7.9 to 16.6 mg/L; and NO₃⁻ from 5.20 to 22.0 mg/L.

The measured stable isotope values for δ¹⁸O ranged from −10.22 to −9.81 ‰, and for δ²H from −71.3 to −69.3 ‰. Preliminary stable isotope results for reservoirs, public taps, and hydrants were plotted against LMWL (Figure 8). It was observed that the summer sampling campaign for reservoirs and public taps differs from the winter sampling (reservoirs and hydrants). This discrepancy can be attributed to different factors such as water temperature variations, fractionations within pipes, different pumping regimes, etc. Regardless of all these influences, it was observed that values are distributed within the elliptical shapes, which represent their origin. Moreover, one reservoir in both sampling campaigns showed the same origin (reservoir R1 and the Bela source). Consistent with the geochemical parameters, there were no outliers observed, confirming no contamination within the water supply network.

3.3. Origin of Water in the Water Supply Network

3.3.1. Origin of Water at Exits at the Bartolovec Wellfield

In the Bartolovec wellfield, mixing modeling was conducted on four of the five exits, as exit Tonimir exclusively receives water from well B-1 (SA). Also, it was observed that modeling with δ¹⁸O as constraints did not show satisfactory matches in most calculations, so these results were not considered. Chloride emerged as the most effective constraint in the majority of MB model calculations, with a generally better fit observed in the IM and FM models compared to MB models (Table 2).

Nevertheless, all three models showed that at exit Doljan, depending on the pumping regime, more than 50% and up to 77% of water is pumped from the DA and distributed towards the southern part of the network. On the other hand, exit Ludbreg primarily supplies the eastern part of the network with DA water, ranging from 92 to 100% (Table 2). At the exit Varaždin, up until April 2022, water was mainly pumped from the DA, between 93 to 99% of the supply. After that, the regime of pumping changed, favoring SA water from 62 to 72%, before reverting back to dominantly DA water (Table 2). At the exit, Novakovec primarily supplies the eastern part of the network with the SA water, ranging from 92 to 100% (Table 2).

3.3.2. Reservoirs, Public Taps and Hydrants

Similar to the exits, the best fitting obtained by MB models is achieved using chloride as constraints, and generally, the IM and FM models display a better fit than MB. After validation, it was observed that IM models yielded the best results, which are discussed further in the text. Despite slightly more difficulties appearing here, given that most of the results showed a good match, it is evident that this modeling approach needs to be further elaborated. The differences exist due to the varied chemical and isotopic composition of many wells, which are alternately integrated into the water supply. It is likely that more accurate results would be obtained by sampling all wells at the water pumping stations rather than just their representatives. In addition, the chemical composition of water changes while standing in reservoirs. During the warm months with higher temperatures, evaporation occurs, which affects the composition of δ¹⁸O, while stagnant water leads to precipitation of carbonates, which changes the EC. Moreover, the amount of free chlorine decreases, and the number of bacteria that can affect nitrification processes increases. Chlorination, albeit to a minor extent, affects chloride concentration, thereby potentially influencing the results. It is possible that each of these processes can affect the modeling result. In addition, by introducing additional indicators of water quality (anions, cations), the accuracy of the models can be increased, as well as by expanding the sampling locations.

A simple Piper diagram showed at a glance that only one reservoir (R1) during both sampling campaigns was filled from only one groundwater source, Bela (Figure 9). However, all other network sampling points demonstrated a mixing of different sources.

It was observed that in the western part of the water supply network, reservoirs predominantly contain a mixture of waters from Bartolovec (exits E1 and E4) and Bela water sources (Figure 10c). Unfortunately, due to the absence of taps in this area, sampling was not conducted. In the northern part of the network, a mixture of waters between Vinokovščak and Bartolovec sources was observed (Figure 10c). Hydrants and taps in closer proximity to the Vinokovščak source had higher ratios, although certain instances of these higher ratios were unexpected. Moreover, differences between summer and winter periods were apparent. The southern and eastern part of the network is mainly filled with water from the Bartolovec source (exits E2 and E3) and B1 (Figure 10a–c). Mixing ratios showed that the main groundwater source is Bartolovec, which has an essential DA role. Examination of pumping amounts during the sampling dates corroborated the modeled values.

4. Conclusions

In a time when water scarcity is an issue not just in developing countries but also in all countries over the globe due to climate change’s impact on groundwater recharge, understanding the water supply network is the key to good water supply management. By using hydrochemical parameters, such as basic anions and physical-chemical properties of water, which are measured in everyday routine in water supplier laboratories, combined with water stable isotopes, extended anion and cation analysis, and geochemical modeling (MB, IM, FM), the water within the water supply network could be comprehensively characterized as shown in this paper. In the Varaždin WSS, all measured parameters are below Maximum Contaminant Level (MCL) values, ensuring the safety of water for human consumption. The differences between water sources were observed based on basic chemical composition. Furthermore, the Piper diagram provided a simplified overview of mixing patterns within the water supply network. Although all three modeling approaches showed satisfactory results, IM results were more reliable than FM and BM. Mixing modeling demonstrated the mixing of all three sources at the western, southwestern, and northern parts of the water supply network. Conversely, the eastern and southeastern parts of the network depend primarily on the DA water from the Bartolovec groundwater source, suggesting that these parts of the network are very vulnerable in the event of an emergency at that source. That should be additionally taken into account during the risk assessment of the water supply system and the establishment of control measures that will further increase the safety of the water supply (enabling the supply of water from other existing sources and increasing the water storage capacity in the threatened area; short-term measure: supplying water by mobile water tanks). Although δ¹⁸O did not provide good mixing results due to fractionation in pipelines, it revealed that water in reservoirs is in equilibrium with air (no evaporation effects), indicating well-sealed and impermeable reservoirs. The results showed that geochemical modeling could be effectively applied in water security management plans. The proposed methodology can be applied globally in a variety of water supply system settings, especially where tremendous water losses are present, e.g., due to leaking pipes, etc.

Author Contributions

Conceptualization, N.N.-H., I.K. (Igor Karlović) and T.M.; methodology T.M.; software, N.N.-H., T.M. and I.K. (Igor Karlović); validation, T.M.; formal analysis, N.N.-H. and T.M.; investigation N.N.-H., I.K. (Igor Karlović) and T.M.; resources N.N.-H. and T.M.; writing—original draft preparation, N.N.-H.; writing—review and editing, T.M., I.K. (Igor Karlović) and I.K. (Ivan Kovač); visualization T.M. and I.K. (Igor Karlović). All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicate.

Informed Consent Statement

Not applicate.

Data Availability Statement

Data are in Supplementary materials.

Acknowledgments

The authors would like to thank the Varaždin Utility Company (VARKOM) and Robert Rutić for their support during sampling campaigns.

Conflicts of Interest

Author Nikolina Novotni-Horčička was employed by the company Varkom d.o.o. The remaining 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 result.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Schematic map of sampled points. See Table 1 for description.

Figure 2. Measured monthly values of (a) EC, (b) Cl−, and (c) NO3− in groundwater of Bartolovec wellfield.

Figure 2. Measured monthly values of (a) EC, (b) Cl−, and (c) NO3− in groundwater of Bartolovec wellfield.

Figure 3. Piper diagram of sampled wells and spring (red: Bela, blue: Zp-7 (SA), purple: Zd-4 (DA), and grey: ZV-4).

Figure 4. Monthly oscillation of δ18O (‰) values in groundwater of the Bartolovec wellfield.

Figure 5. Measured monthly values of EC, Cl− and NO3− in groundwater of Vinokovščak wellfield.

Figure 6. Distribution of groundwater δ18O (‰) and δ2H (‰) values around the local meteoric water line (LMWL) (according to [11]).

Figure 7. Measured monthly values of EC, Cl− and NO3− in groundwater of Bela spring.

Figure 8. The relationship between δ2H and δ18O in drinking water supply network, plotted alongside the LMWL (according to Marković et al., 2020 [11]). The range of measured stable isotopes for groundwater sources is illustrated in the form of elliptical shapes, with each color representing a different source: blue—Bela; red—Vinokovščak; purple—DA Bartolovec; grey—SA Bartolovec.

Figure 9. Piper diagram of sampled waters from the network (brown color—reservoirs; green color—public taps; blue color—hydrants).

Figure 10. Mixing ratios in (a) public taps, (b) hydrants, and (c) reservoirs.

Table 1

Sampling points coordinates and technical data.

Mark on Map (Figure 1)	Code	Description	Longitude (E) (HTRS96/TM)	Latitude (N) (HTRS96/TM)	Elevation(m a.s.l.)	Well Depth (m)	Filter (m)
1	H1	Hydrant AGM	487,061	5,130,938	172.24
2	H2	Hydrant Perivoj M.Culeja	486,712	5,129,729	171.72
3	H3	Hydrant, S. Vukovića	487,617	5,128,649	170.00
4	H4	Hydrant Drava	488,490	5,131,493	170.50
5	H5	Hydrant Sračinec	483,780	5,131,782	175.12
6	H6	Hydrant Ludbreg	509,239	5,123,187	153.50
7	H7	Hydrant Varaždinske Toplice	494,099	5,118,929	189.72
8	PT1	Public tap AGM	487,061	5,130,938	172.24
9	PT2	Public tap Perivoj M.Culeja	486,712	5,129,729	171.72
10	PT3	Public tap S. Vukovića	487,609	5,128,958	171.54
11	PT4	Public tap Drava	488,490	5,131,493	170.50
12	PT5	Public tap Sračinec	483,780	5,131,782	153.50
13	PT6	Public tap Ludbreg	509,239	5,123,187	175.12
14	PT7	Public tap Varaždinske Toplice	494,156	5,118,958	192.23
15	R1	Reservoir Briška	483,631	5,121,191	331.48
16	R2	Reservoir Bolfan	513,429	5,118,567	258.52
17	R3	Reservoir Breznički Hum	483,103	5,107,850	240.00
18	R4	Reservoir Doljan	487,150	5,122,584	223.13
19	R5	Reservoir Donja Voća	470,694	5,130,305	275.30
20	R6	Reservoir Lepoglava	465,054	5,119,159	285.98
21	R7	Reservoir Ludbreg	508,931	5,122,285	207.58
22	R8	Reservoir Lužan	489,003	5,119,909	305.41
23	R9	Reservoir Pintarići	465,842	5,128,636	432.18
24	R10	Reservoir Šćepanje	482,649	5,109,059	336.34
25	R11	Reservoir Tonimir	492,999	5,120,354	258.81
26	R12	Reservoir Topličica	483,721	5,114,068	358.96
27	R13	Reservoir Vinica	472,327	5,132,738	259.86
28	B1	Well B-1, Bartolovec well field (BWF), SA	494,829	5,127,376	161.56	32	15–25
29	Zd-4	Well Zd-4, BWF, DA	494,858	5,127,550	162.50	105	62–102
30	Zp-7	Well Zp-7, BWF, SA	495,036	5,127,384	161.35	45	15–42
31	Zp-9	Well Zp-9, BWF, SA	494,671	5,127,620	162.14	45	15–42
32	E1	Exit Varaždin, BWF	494,880	5,127,604	162.13
33	E2	Exit Ludbreg, BWF	494,880	5,127,604	162.13
34	E3	Exit Novakovec, BWF	494,978	5,127,382	162.13
35	E4	Exit Doljan, BWF	494,832	5,127,378	162.13
36	ZV4	Well ZV4, Vinokovščak well field	485,842	5,133,689	174.80	45	10–40
37	Bela	Bela spring	480,891	5,117,696	218.54

Table 2

Results of MB, FM, and IM modeling evaluated against measured values at all exits. Red color—not good matching; yellow color—poor matching; green color—excelent matching.

exit Doljan	MB model 1		MB model 2		MB model 3		FM model		IM model		Measured	Calculated MB1	Calculated MB2	Calculated MB3	Calculated FM	Calculated IM	Calculated MB 1	Calculated MB 2	Calculated MB 3	Calculated FM	Calculated IM
	DA	SA	DA	SA	DA	SA	DA	SA	DA	SA	SO42− (mg/L)	SO42− (mg/L)	SO42− (mg/L)	SO42− (mg/L)	SO42− (mg/L)	SO42− (mg/L)	deviation	deviation	deviation	deviation	deviation
Dec-21	78%	22%	71%	29%	73%	27%	71%	29%	72%	28%	12.4	11.6	13.3	12.8	13.3	13.0	0.17	1.04	0.52	1.05	0.70
Jan-22	80%	20%	76%	24%	69%	31%	76%	24%	77%	23%	15.2	14.4	15.5	17.3	15.5	15.2	0.16	0.35	2.55	0.38	0.06
Feb-22	77%	23%	73%	27%	70%	30%	73%	27%	74%	26%	16.2	15.6	16.7	17.4	16.7	16.4	0.13	0.61	1.45	0.63	0.28
Mar-22	77%	23%	75%	25%	74%	26%	75%	25%	75%	25%	13.1	13.1	13.6	14.0	13.6	13.7	0.00	0.63	1.04	0.61	0.73
Apr-22	78%	22%	76%	24%	79%	21%	76%	24%	76%	24%	13.5	13.1	13.7	12.9	13.7	13.7	0.08	0.30	0.13	0.30	0.20
May-22	67%	33%	58%	42%	85%	15%	58%	42%	59%	41%	15.7	12.8	14.6	9.1	14.6	14.4	0.61	0.23	1.36	0.24	0.28
Jun-22	70%	30%	59%	41%	85%	15%	59%	41%	61%	39%	16.0	12.7	14.8	9.7	14.8	14.5	0.69	0.25	1.31	0.25	0.32
Jul-22	69%	31%	60%	40%	79%	21%	60%	40%	60%	40%	16.9	13.9	16.1	11.7	16.2	16.2	0.61	0.16	1.07	0.15	0.14
Aug-22	80%	20%	67%	33%	79%	21%	69%	31%	71%	29%	14.1	11.2	14.4	11.5	13.9	13.5	0.60	0.36	0.55	0.04	0.13
Sep-22	79%	21%	67%	33%	78%	22%	70%	30%	71%	29%	14.0	11.7	14.7	12.0	13.9	13.7	0.47	0.84	0.41	0.03	0.07
Oct-22	64%	36%	61%	39%	74%	26%	55%	45%	58%	42%	17.5	15.0	15.7	12.6	17.2	16.5	0.51	0.37	1.02	0.07	0.20
Nov-22	65%	35%	59%	41%	75%	25%	56%	44%	59%	41%	16.5	14.5	15.8	12.1	16.6	15.9	0.42	0.14	0.92	0.13	0.11
Dec-22	75%	25%	62%	38%	62%	38%	64%	36%	65%	35%	14.7	12.2	15.5	15.3	15.0	14.7	0.51	0.93	0.77	0.37	0.01
																std	0.74	0.30	0.61	0.29	0.23
exit Ludbreg	MB model 1		MB model 2		MB model 3		FM model		IM model		Measured	Calculated MB1	Calculated MB2	Calculated MB3	Calculated FM	Calculated IM	Calculated MB 1	Calculated MB 2	Calculated MB 3	Calculated FM	Calculated IM
	DA	SA	DA	SA	DA	SA	DA	SA	DA	SA	SO42− (mg/L)	SO42− (mg/L)	SO42− (mg/L)	SO42− (mg/L)	SO42− (mg/L)	SO42− (mg/L)	deviation	deviation	deviation	deviation	deviation
Dec-21	99%	1%	93%	7%	99%	1%	93%	7%	99%	1%	7.2	6.3	7.8	6.4	7.8	6.4	0.18	0.75	0.17	0.75	0.16
Jan-22	100%	0%	99%	1%	98%	2%	99%	1%	98%	2%	9.1	9.2	9.5	9.8	9.4	9.8	0.11	0.38	0.76	0.36	0.74
Feb-22	100%	0%	97%	3%	99%	1%	97%	3%	99%	1%	9.6	9.6	10.4	9.9	10.4	9.9	0.02	0.93	0.29	0.94	0.34
Mar-22	99%	1%	98%	2%	105%	−5%	98%	2%	100%	0%	7.1	7.5	7.6	6.0	7.6	7.2	0.46	0.62	0.23	0.61	0.11
Apr-22	100%	0%	98%	2%	100%	0%	98%	2%	100%	0%	7.2	7.3	7.7	7.2	7.7	7.3	0.10	0.58	0.01	0.61	0.10
May-22	100%	0%	99%	1%	100%	0%	100%	1%	100%	0%	5.4	6.0	6.1	5.9	6.1	6.0	0.71	0.86	0.60	0.84	0.71
Jun-22	71%	29%	98%	2%	100%	0%	98%	2%	100%	0%	5.3	12.4	6.9	6.5	7.0	6.5	8.52	1.98	1.40	2.01	1.47
Jul-22	71%	29%	100%	0%	100%	0%	100%	0%	100%	0%	6.7	13.7	6.5	6.4	6.5	6.5	8.45	0.03	0.06	0.03	0.04
Aug-22	100%	0%	100%	0%	99%	1%	100%	0%	99%	1%	6.7	6.5	6.5	6.6	6.5	6.7	0.04	0.04	0.01	0.04	0.05
Sep-22	100%	0%	99%	1%	99%	1%	99%	1%	99%	1%	6.3	6.8	7.1	7.0	7.1	7.0	0.60	0.98	0.79	0.94	0.89
Oct-22	100%	0%	95%	5%	99%	1%	94%	6%	97%	3%	6.5	6.2	7.4	6.5	7.6	6.9	0.06	1.11	0.04	1.34	0.53
Nov-22	100%	0%	97%	3%	100%	0%	97%	3%	100%	0%	5.0	5.9	6.6	5.9	6.6	5.9	1.09	1.91	1.09	1.97	1.09
Dec-22	100%	0%	99%	1%	92%	8%	99%	1%	94%	6%	6.3	6.2	6.5	8.2	6.5	7.7	0.02	0.26	2.33	0.24	1.63
																std	3.09	0.46	0.69	0.44	0.50
exit Varaždin	MB model 1		MB model 2		MB model 3		FM model		IM model		Measured	Calculated MB1	Calculated MB2	Calculated MB3	Calculated FM	Calculated IM	Calculated MB 1	Calculated MB 2	Calculated MB 3	Calculated FM	Calculated IM
	DA	SA	DA	SA	DA	SA	DA	SA	DA	SA	SO42− (mg/L)	SO42− (mg/L)	SO42− (mg/L)	SO42− (mg/L)	SO42− (mg/L)	SO42− (mg/L)	deviation	deviation	deviation	deviation	deviation
Dec-21	100%	0%	100%	0%	71%	29%	99%	1%	99%	1%	8.2	6.2	6.1	13.3	6.4	6.4	0.41	0.44	6.19	0.36	0.36
Jan-22	100%	0%	94%	6%	70%	30%	94%	6%	95%	5%	10.7	9.2	10.8	17.2	10.8	10.5	0.30	0.11	7.79	0.09	0.03
Feb-22	100%	0%	93%	7%	62%	38%	93%	7%	93%	7%	10.7	9.6	11.5	19.6	11.5	11.5	0.22	0.99	10.79	0.99	0.92
Mar-22	99%	1%	94%	6%	62%	38%	94%	6%	94%	6%	8.5	7.4	8.9	17.1	8.9	8.8	0.23	0.48	10.34	0.48	0.34
Apr-22	100%	0%	93%	7%	49%	51%	93%	7%	93%	7%	8.7	7.3	9.0	20.7	9.1	9.2	0.29	0.37	14.47	0.39	0.52
May-22	73%	27%	31%	69%	92%	8%	31%	69%	31%	69%	21.9	11.6	20.1	7.7	20.1	20.1	2.14	0.37	2.94	0.37	0.38
Jun-22	100%	0%	31%	69%	90%	10%	31%	69%	32%	68%	22.8	6.5	20.7	8.6	20.7	20.4	3.37	0.44	2.95	0.44	0.50
Jul-22	100%	0%	31%	69%	87%	13%	31%	69%	33%	67%	24.8	6.5	23.4	9.7	23.4	22.8	3.80	0.30	3.12	0.30	0.41
Aug-22	87%	13%	63%	37%	78%	22%	66%	34%	63%	37%	15.6	9.6	15.3	11.9	14.7	15.4	1.23	0.06	0.77	0.18	0.05
Sep-22	85%	15%	62%	38%	75%	25%	66%	34%	62%	38%	15.6	10.3	15.8	12.8	14.8	15.8	1.11	0.22	0.58	0.16	0.25
Oct-22	65%	35%	27%	73%	82%	18%	16%	84%	28%	72%	25.5	14.8	24.1	10.6	26.9	23.9	2.22	0.28	3.10	1.68	0.33
Nov-22	84%	16%	67%	33%	81%	19%	64%	36%	67%	33%	15.0	9.8	14.1	10.6	14.7	14.0	1.07	0.19	0.90	0.06	0.21
Dec-22	84%	16%	61%	39%	63%	37%	63%	37%	65%	35%	14.6	10.1	15.6	15.2	15.2	14.7	0.93	1.24	0.69	0.67	0.08
																std	1.20	0.34	4.56	0.44	0.24
exit Novakovec	MB model 1		MB model 2		MB model 3		FM model		IM model		Measured	Calculated MB1	Calculated MB2	Calculated MB3	Calculated FM	Calculated IM	Calculated MB 1	Calculated MB 2	Calculated MB 3	Calculated FM	Calculated IM
	DA	SA	DA	SA	DA	SA	DA	SA	DA	SA	SO42− (mg/L)	SO42− (mg/L)	SO42− (mg/L)	SO42− (mg/L)	SO42− (mg/L)	SO42− (mg/L)	deviation	deviation	deviation	deviation	deviation
Dec-21	1%	99%	7%	93%	0%	100%	65%	35%	31%	69%	29.8	30.2	28.6	30.4	34.5	22.9	0.42	0.24	0.72	5.62	1.43
Jan-22	1%	99%	5%	95%	0%	100%	39%	61%	16%	84%	35.2	35.1	34.1	35.4	39.0	31.2	0.02	0.24	0.24	4.58	0.83
Feb-22	0%	100%	6%	94%	0%	100%	53%	47%	11%	89%	35.0	35.8	34.4	35.8	40.9	32.9	0.97	0.13	0.97	7.14	0.43
Mar-22	0%	100%	9%	91%	0%	100%	100%	0%	6%	94%	33.1	33.2	30.8	33.2	40.4	31.6	0.12	0.48	0.12	8.82	0.30
Apr-22	0%	100%	4%	96%	20%	80%	57%	44%	4%	96%	33.6	33.8	32.7	28.5	37.9	32.7	0.24	0.18	1.06	5.23	0.18
May-22	0%	100%	3%	97%	0%	100%	64%	36%	33%	67%	28.1	26.4	25.8	26.4	30.2	19.7	0.35	0.48	0.35	2.55	1.75
Jun-22	0%	100%	6%	94%	0%	100%	74%	26%	21%	79%	29.5	26.9	25.7	26.9	31.7	22.6	0.54	0.78	0.54	2.67	1.42
Jul-22	2%	98%	6%	94%	0%	100%	6%	94%	9%	91%	30.9	30.5	29.4	30.9	31.3	28.7	0.08	0.32	0.00	0.48	0.46
Aug-22	8%	92%	12%	88%	0%	100%	12%	88%	10%	90%	30.6	28.5	27.6	30.5	31.3	28.1	0.44	0.62	0.02	0.82	0.52
Sep-22	7%	93%	10%	90%	0%	100%	83%	17%	8%	92%	30.5	28.8	28.2	30.5	36.1	28.6	0.35	0.48	0.00	6.78	0.39
Oct-22	13%	87%	11%	89%	0%	100%	0%	100%	10%	90%	30.4	27.6	28.1	30.8	30.8	28.3	0.58	0.48	0.48	0.48	0.43
Dec-22	12%	88%	9%	91%	0%	100%	2%	98%	9%	91%	30.2	27.4	28.3	30.4	30.5	28.2	0.58	0.40	0.24	0.39	0.41
																std	0.26	0.19	0.37	2.96	0.50

Supplementary Materials

The following supporting information can be downloaded at https://www.mdpi.com/article/10.3390/su16219558/s1, Table S1: Results of hydrochemical analysis of sampled groundwater sources; Table S2: Results of hydrochemical analysis of public taps, hydrant, and reservoirs.

References

1. WHO Fact Sheets, Drinking Water. Available online: https://www.who.int/news-room/fact-sheets/detail/drinking-water (accessed on 15 February 2023).

2. Gleick, P.H. Global freshwater resources: Soft-path solutions for the 21st century. Science; 2003; 302, pp. 1524-1528. [DOI: https://dx.doi.org/10.1126/science.1089967] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/14645837]

3. Tipple, B.J.; Jameel, Y.; Chau, T.H.; Mancuso, C.J.; Bowen, G.J.; Dufour, A.; Chesson, L.A.; Ehleringer, J.R. Stable hydrogen and oxygen isotopes of tap water reveal structure of the San Francisco Bay Area’s water system and adjustments during a major drought. Water Res.; 2017; 119, pp. 212-224. [DOI: https://dx.doi.org/10.1016/j.watres.2017.04.022] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28463769]

4. Directive (EU) 2020/2184 on the Quality of Water Intended for Human Consumption. Available online: https://www.fao.org/faolex/results/details/en/c/LEX-FAOC201243/ (accessed on 29 June 2024).

5. Puust, R.; Kapelan, Z.; Savic, D.A.; Koppel, T. A review of methods for leakage management in pipe networks. Urban Water J.; 2013; 7, pp. 25-45. [DOI: https://dx.doi.org/10.1080/15730621003610878]

6. Mutikanga, H.E.; Sharma, S.K.; Vairavamoorthy, K. Methods and tools for managing losses in water distribution systems. J. Water Resour. Plan. Manag.; 2013; 139, pp. 166-174. [DOI: https://dx.doi.org/10.1061/(ASCE)WR.1943-5452.0000245]

7. Waldrip, S.H.; Niven, R.K.; Abel, M.; Schlegel, M. Maximum Entropy Analysis of Hydraulic Pipe Flow Networks. J. Hydraul. Eng.; 2016; 142, 04016028. [DOI: https://dx.doi.org/10.1061/(ASCE)HY.1943-7900.0001126]

8. Froehlich, K.; Gibson, J.J.; Aggarwal, P.K. Deuterium Excess in Precipitation and Its Climatological Significance. No. IAEA-CSP--13/P 2002; Available online: https://inis.iaea.org/search/search.aspx?orig_q=RN:34017972 (accessed on 29 June 2024).

9. Dotsika, E.; Lykoudis, S.; Poutoukis, D. Spatial distribution of the isotopic composition of precipitation and spring water in Greece. Global Planet Chang.; 2010; 71, pp. 141-149. [DOI: https://dx.doi.org/10.1016/j.gloplacha.2009.10.007]

10. Gross, E.; Andrews, S.; Bergamaschi, B.; Downing, B.; Holleman, R.; Burdick, S.; Durand, D. The Use of Stable Isotope-Based Water Age to Evaluate a Hydrodynamic Model. Water; 2019; 11, 2207. [DOI: https://dx.doi.org/10.3390/w11112207]

11. Marković, T.; Karlović, I.; Perčec Tadić, M.; Larva, O. Application of Stable Water Isotopes to Improve Conceptual Model of Alluvial Aquifer in the Varaždin Area. Water; 2020; 12, 379. [DOI: https://dx.doi.org/10.3390/w12020379]

12. Vreča, P.; Kern, Z. Use of water isotopes in hydrological processes. Water; 2020; 12, 2227. [DOI: https://dx.doi.org/10.3390/w12082227]

13. Jameel, Y.; Brewer, S.; Good, S.P.; Tipple, B.J.; Ehleringer, J.R.; Bowen, G.J. Tap water isotope ratios reflect urban water system structure and dynamics across a semiarid metropolitan area. Water Resour. Res.; 2016; 52, pp. 5891-5910. [DOI: https://dx.doi.org/10.1002/2016WR019104]

14. Yusuf, J.; Simon, B.; Richard, P.F.; Brett, J.T.; Shazelle, T.; Gabriel, J.B. Isotopic reconnaissance of urban WSS dynamics. Hydrol. Earth Syst. Sci.; 2018; 22, pp. 6109-6125. [DOI: https://dx.doi.org/10.5194/hess-22-6109-2018]

15. Grimmeisen, F.; Lehmann, M.F.; Liesch, T.; Goeppert, N.; Klinger, J.; Zopfi, J.; Goldscheider, N. Isotopic constraints on water source mixing, network leakage and contamination in an urban groundwater system. Sci. Total Environ.; 2017; 583, pp. 202-213. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2017.01.054] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28104328]

16. Nagode, K.; Kanduč, T.; Zuliani, T.; Bračič Železnik, B.; Jamnik, B.; Vreča, P. Daily Fluctuations in the Isotope and Elemental Composition of Tap Water in Ljubljana, Slovenia. Water; 2021; 13, 1451. [DOI: https://dx.doi.org/10.3390/w13111451]

17. Kurajica, L.; Ujević Bošnjak, M.; Kinsela, A.S.; Štiglić, J.; Waite, T.D.; Capak, K.; Pavlić, Z. Effects of changing supply water quality on drinking water distribution networks: Changes in NOM optical properties, disinfection byproduct formation, and Mn deposition and release. Sci. Total Environ.; 2021; 762, 144159. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2020.144159] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33360458]

18. Kurajica, L.; Ujević Bošnjak, M.; Kinsela, A.S.; Bieroza, M.; Štiglić, J.; Waite, T.D.; Capak, K.; Romić, Ž. Mixing of arsenic-rich groundwater and surface water in drinking water distribution systems: Implications for contaminants, disinfection byproducts and organic component. Chemosphere; 2022; 292, 133406. [DOI: https://dx.doi.org/10.1016/j.chemosphere.2021.133406] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34958791]

19. Horvatinčić, N.; Kapelj, S.; Sironić, A.; Krajcar Bronić, I.; Kapelj, J.; Marković, T. Investigation of water resources and water protection in the karst area of Croatia using isotopic geochemical analyses. Proceedings of the Advances in Isotope Hydrology and its Role in Sustainable Water Resources Management, Book of Abstracts IAEA-CN-151; Wiena, Austria, 21–25 May 2007; Volume 2, pp. 295-304.

20. Buškulić, P.; Parlov, J.; Kovač, Z.; Brenko, T.; Pejić, M. Determination of Nitrate Migration and Distribution through Eutric Cambisols in an Area without Anthropogenic Sources of Nitrate (Velika Gorica Well Field, Croatia). Sustainability; 2023; 15, 16529. [DOI: https://dx.doi.org/10.3390/su152316529]

21. Cukrov, N.; Tepić, N.; Omanović, D.; Lojen, S.; Bura-Nakić, E.; Vojvodić, V.; Pižeta, I. Qualitative interpretation of physico-chemical and isotopic parameters in the Krka River (Croatia) assessed by multivariate statistical analysis. JIEC; 2012; 92, pp. 1187-1199. [DOI: https://dx.doi.org/10.1080/03067319.2010.550003]

22. Ozyurt, N.N.; Lutz, H.O.; Hunjak, T.; Mance, D.; Roller-Lutz, Z. Characterization of the Gacka River basin karst aquifer (Croatia): Hydrochemistry, stable isotopes and tritium-based mean residence times. Sci. Total Environ.; 2014; 487, pp. 245-254. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2014.04.018] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24784749]

23. Mance, D.; Radišić, M.; Lenac, D.; Rubinić, J. Hydrological Behavior of Karst Systems Identified by Statistical Analyses of Stable Isotope Monitoring Results. Hydrology; 2022; 9, 82. [DOI: https://dx.doi.org/10.3390/hydrology9050082]

24. Parlov, J.; Kovač, Z.; Barešić, J. The study of the interactions between Sava River and Zagreb aquifer system (Croatia) using water stable isotopes. Proceedings of the E3S Web of Conferences; Tomsk, Russia, 21–26 July 2019; Volume 98, 12017.

25. Krajcar Bronić, I.; Barešić, J.; Sironić, A.; Lovrenčić Mikelić, I.; Borković, D.; Horvatinčić, N.; Kovač, Z. Isotope Composition of Precipitation, Groundwater, and Surface and Lake Waters from the Plitvice Lakes, Croatia. Water; 2020; 12, 2414. [DOI: https://dx.doi.org/10.3390/w12092414]

26. Mance, D.; Lenac, D.; Rubinić, J. Isotope studies of karst springs included in the water supply system of the City of Rijeka (Croatia). Southeast. Eur. Med. J.; 2017; 1, pp. 46-54.

27. Urumović, K.; Hlevnjak, B.; Prelogović, E.; Mayer, D. Hydrogeological conditions of Varaždin aquifer. Geol. Vjesn.; 1990; 43, pp. 149-158.

28. Crnogaj, S.; Mraz, V.; Kruk, B.; Jurić, A. Study on the reserve technical building stone dolomite in mine “Belski Dol Gornji” near Podrute; Croatian Geological Survey: Zagreb, Croatia, 1997; (In Croatian)

29. Mandel, S.; Shiftan, Z.L. Groundwater Resources Investigation and Development; Academic Press: New York, NY, USA, 1981.

30. Domenico, P.A.; Schwartz, F.W. Physical and Chemical Hydrogeology; John Willey and Sons: New York, NY, USA, 1990.

31. Picarro Technology. Cavity Ring down Spectroscopy. Available online: http://www.picarro.com/technology/cavity_ring_down_spectroscopy (accessed on 24 September 2023).

32. Clark, I.D.; Fritz, P. Environmental Isotopes in Hydrogeology; CRC Press: Boca Raton, FL, USA, 2013.

33. Craig, H. Isotopic Variations in Meteoric Waters. Science; 1961; 133, pp. 1702-1703. [DOI: https://dx.doi.org/10.1126/science.133.3465.1702] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17814749]

34. Vázquez-Suñé, E.; Carrera, J.; Tubau, I.; Sánchez-Vila, X.; Soler, A. An approach to identify urban groundwater recharge. Hydrol. Earth Syst. Sci.; 2010; 14, pp. 2085-2097. [DOI: https://dx.doi.org/10.5194/hess-14-2085-2010]

35. Eberts, S.M.; Thomas, M.A.; Jagucki, M.L. Factors Affecting Public-Supply Well Vulnerability to Contamination: Understanding Observed Water Quality and Anticipating Future Water Quality; US Department of the Interior: Washington, DC, USA, US Geological Survey: Reston, VA, USA, 2013; 132p.

36. Hepburn, E.; Cendón, D.I.; Bekele, D.; Currell, M. Environmental isotopes as indicators of groundwater recharge, residence times and salinity in a coastal urban redevelopment precinct in Australia. Hydrogeol. J.; 2020; 28, pp. 503-520. [DOI: https://dx.doi.org/10.1007/s10040-019-02077-x]

37. Parkhurst, D.L.; Appelo, C.A.J. Description of Input and Examples for PHREEQC Version 3—A Computer Program for Speciation, Batch-Reaction, One-Dimensional Transport, and Inverse Geochemical Calculations; U.S. Geological Survey Techniques and Methods, Book 6 USGS: Denver, CO, USA, 2013; Chapter A43 497.Available online: https://pubs.usgs.gov/tm/06/a43/ (accessed on 23 September 2023).

38. Plummer, L.N.; Prestemon, E.C.; Parkhurst, D.L. NETPATH: An Interactive Code (NETPATH) for Modelling NET Geochemical Reactions along a Flow PATH, Version 2.0; Water-Resources Investigations Report 94-4169 USGS: Reston, VA, USA, 1994.

39. El-Kadi, A.I.; Plummer, L.N.; Aggarwal, P. NETPATH-WIN: An interactive user version of the mass-balance model, NETPATH. Ground Water; 2010; 49, pp. 593-599. [DOI: https://dx.doi.org/10.1111/j.1745-6584.2010.00779.x]

40. Kovač, I.; Kovačev-Marinčić, B.; Novotni-Horčička, N.; Mesec, J.; Vugrinec, J. Komparativna analiza koncentracije nitrata u gornjem i donjem sloju Varaždinskog vodonosnika. Radovi Zavoda za znanstveni rad Varaždin; 2017; 28, pp. 41-57. (In Croatian) [DOI: https://dx.doi.org/10.21857/9e31lh4zem]

41. Karlović, I.; Marković, T.; Šparica Miko, M.; Maldini, K. Geochemical Characteristics of Alluvial Aquifer in the Varaždin Region. Water; 2021; 13, 1508. [DOI: https://dx.doi.org/10.3390/w13111508]

42. Babić, Ž.; Čakarun, I.; Sokač, A.; Mraz, V. On geological features of quaternary sediments of Drava basin on Croatian territory. Geol Vjesn; 1978; 30, pp. 43-61.

43. Urumović, K. On quaternary aquifer complex in Varaždin area. Geol. Vjesn.; 1971; 24, pp. 183-191.

44. Marković, T.; Larva, O. Hydrogeological characterization of carbonate aquifer system in the northwestern part of Croatia, Abstracts Book. Proceedings of the IAH 40th International Congress; Perth, Australia, 15 September–20 September 2013; 179.

45. Zaninović, K.; Gajić-Čapka, M.; Perčec Tadić, M.; Vučetić, M.; Milković, J.; Bajić, A.; Cindrić, K.; Cvitan, L.; Katušin, Z.; Kaučić, D. et al. Klimatski atlas Hrvatske/Climate atlas of Croatia 1961–1990, 1971–2000; Meteorological and Hydrological Service of Croatia: Zagreb, Croatia, 2008; 200.