Overview of the Chemical and Isotopic Investigations of the Mareza Springs and the Zeta River in Montenegro

Abstract

The Mareza karst aquifer is the most important drinking water resource for the water supply system of the City of Podgorica, the capital of Montenegro. This study presents the first assessment for the determination of the Mareza catchment area. Water chemistry and stable isotopic composition (δ¹⁸O and δ²H) of monthly precipitation samples (as inputs) are presented, in order to determine the Local Meteoric Water Line (LMWL) for the study area, and to analyze the behavior of the karst spring Mareza (as output) and the Zeta River water. The possible impact of the river on the Mareza springs was also investigated. Stable isotope compositions were used to analyze the origin of the four springs of the Mareza aquifer. Seasonal variations of δ¹⁸O and δ²H values and deuterium excess (d excess) changes in precipitation are explained by the mixing of air masses, such that a Mediterranean source prevails in the winter period, while in the summer period, the area is rather under the influence of air mass originating from the Atlantic Ocean. All spring water samples have lower δ values than the local precipitation and they plot above the LMWL, which may indicate recharge at a higher altitude in the distant mountainous area. The d excess values of all water samples (higher than 10‰) indicate the prevalence of the Mediterranean as a moisture source. Based on the analysis of the seasonal variations of δ¹⁸O and δ²H in precipitation and the Mareza spring, it has been estimated that the groundwater mean transit time (MTT) is 92–129 days, and that the young water fraction (F_yw) amounts to 40.9%–53.3%. These values are typical for the strong karstic springs of highly karstified terrains.

1. Introduction

Karst aquifers are of utmost importance for the drinking water supply in Montenegro, where more than two-thirds of the territory are covered by carbonate formations (limestone, dolostone, and dolomite) belonging to south-eastern Dinaric Alps, and more than 90% of the population depends upon drinking water supplied by groundwater. Rapid development of tourism in the last two decades has increased the need for drinking water, in particular in the coastal areas, during the summer, where the population during the peak season has more than tripled [1].

Compared to other groundwater systems, karstic aquifers are intrinsically vulnerable, since karstified carbonate terrains exhibit strong fissuring and spatially highly variable multiphase porosity, which enables diffuse and concentrated infiltration of water from the surface. Karstic groundwater flow is therefore different from that in other types of aquifers [2,3], namely a combination of slow fracture and quick conduit flows, where the water residence time in the matrix can be several orders of magnitude longer than that in the conduits [4]. Infiltration rates and underground water flow are controlled by heterogeneous permeability and hydrological conditions. The abundance and quality of karstic water resources is thus highly susceptible to climate [5] and land use change [6].

Protection and sustainable management of karstic groundwater require sound evaluation of recharge based on adequate delineation of recharge areas, which is regularly supported by the use of a variety of natural or artificial geochemical and isotopic tracers [7,8,9]. The Mareza aquifer, which provides drinking water for more than 60% of the municipal area of Podgorica (120,000 inhabitants), has been the subject of many scientific investigations [10,11,12,13,14,15], but until now, isotopic tracers have not been used. The origin of the Mareza spring and the nearby Zeta River water was interpreted on the basis of hydrogeological characteristics and structure of the terrain. The results of a previous study of the Mareza were inconclusive, since it could not be identified whether it is a uniform or seasonally fragmented karst aquifer or a group of aquifers which do not communicate with each other [14]. A general direction of the water flow from north to south was identified, but the recharge area has not yet been determined.

Stable O and H isotopes are excellent natural tracers of sources of water because of the spatial and temporal variability of the O and H isotope ratios in precipitation [16]. They have greatly contributed to the understanding of hydrological processes [17,18], particularly that of groundwater recharge [19,20,21], groundwater mixing [22,23,24,25], and the interactions between groundwater and surface water [26,27,28,29]. On the other hand, groundwater chemistry is controlled greatly by interaction of water with the aquifer through which it flows. In any area, groundwater exhibits an isotopic and chemical fingerprint, obtained during processes such as soil/rock-water interaction during recharge and groundwater flow, prolonged storage in the aquifer, dissolution and precipitation of minerals, element speciation, and so forth [30]. The chemistry of karstic ground and spring water and the temporal and spatial variations of major ion concentrations are important in understanding chemical and physical processes that are influenced by geological setting, climate, and anthropogenic activities [31]. These processes can also be studied by determining the spatial and temporal variability of the isotopic composition of hydrogen, oxygen, and carbon in water [32]. Stable isotopes of hydrogen and oxygen are useful tools to indicate the potential water inputs to the system and to determine the mixing of waters with different origins within the system [33,34].

In the present paper, we report the first results of water isotope analyses of precipitation and the Zeta River and Mareza spring water. Based on the results, the origin and mean transit time of spring and river water were determined, providing essential information for groundwater management in the area.

2. Materials and Methods

2.1. Study Area

The Mareza karst aquifer is located in the central part of Montenegro (Figure 1). It belongs to the Adriatic hydrological basin—more precisely, the Skadar Lake basin—and lays about 50 km away from the Adriatic Sea and 20 km away from the lake Skadar. The Mareza springs are located at elevations from 26 to 32.7 m asl (above sea level). The Zeta River, with a total length of 85 km, flows a few kilometers northeast of the Mareza springs.

The study area includes the lower course of the Zeta River (the valley is known as the Zeta River Valley or the Bjelopavlici Valley), its permanent or ephemeral tributaries, and springs. The valley is the only traffic corridor between the two biggest Montenegrin towns of Nikšić and Podgorica. The Mareza aquifer and the Zeta River are separated by the Velje hill, which is built of permeable limestone (Figure 1 and Figure 2). The area around the Mareza springs consists of highly permeable karstified Upper Cretaceous limestone, which lays on clayey sediments and clays [35]. Generally, the recharge area of the Skadar Lake karst aquifer is spread over a wide area in the karst plateaus, while the discharge area extends along the karst depressions and deep canyons [36].

According to the Institute of Hydrometeorology and Seismology of Montenegro’s (IHMS) meteorological data, the mean monthly air temperature during the sampling period (from February 2017 to March 2018) in Podgorica was between 7.3 °C in December and 30 °C in July and August (Figure 3). No data for the town of Danilovgrad were available. The lowest daily air temperature was recorded in December (−3.3 °C) and the highest in July and August (up to 43 °C). The annual average temperature (16.1 °C) was slightly lower than the long-term (2007–2018) average annual temperature (16.8 °C).

During the sampling period, the highest total monthly precipitation in Podgorica occurred in the winter, reaching up to 460 mm (Figure 3). Extreme rainfall events occurred in November 2017 (13 November, 114.5 mm; 30 November, 75.3 mm) and on 3 February 2018 (88 mm), which caused increased turbidity of the water in both the springs and the river. The situation was similar in the town of Danilovgrad. The lowest amount of precipitation occurred in summer (June, July, and August; Figure 3). The total annual precipitation during the sampling period in the study area was 2153 mm for Danilovgrad and 1905 mm for Podgorica.

For comparison, recorded data of mean annual precipitation across a long-term period (2007–2017) in the study area (1773 mm) and data in the study period (1905 mm) showed no significant changes (Supplementary Materials, Figure S1).

2.2. Sampling Locations

Precipitation was collected in rain gauge HDPE (high-density polyethylene) collectors as cumulative monthly samples, positioned at two locations: P1 (near the Zeta River) and P2 (between the capped Mareza springs). Both sites are very close to the hydrometeorological monitoring stations of Danilovgrad (P1) and Podgorica (P2) managed by the IHMS (

Table 1. Geographical latitude, longitude, and altitude of the precipitation (P1 and P2) sampling points and the meteorological stations.

Location	Latitude (°)	Longitude (°)	Altitude (m asl)
P1	42.36	19.04	38
P2	42.45	19.20	32
Danilovgrad meteorological station	42.55	19.1	53
Podgorica meteorological station	42.43	19.28	49

and Figure 1).

Water samples for chemical and isotope analyses were taken at two locations from the Zeta River (upstream, RU; downstream, RD; Figure 1) and from four springs of the Mareza aquifer, of which three are capped (S1, S2, and S3) and one is an open spring (S4). Sampling point RD is located approximately 20 km downstream from RU.

2.3. Sampling

During the one-year sampling period, 74 spring water samples, 37 river water samples, and 21 samples of precipitation were collected and analyzed.

Precipitation was collected in 5-L HDPE collectors. To prevent evaporation, even under very hot summer conditions, a layer of paraffin oil was added [37]. Before sample bottles were filled, the paraffin oil was removed completely using a glass separation funnel. Untreated precipitation samples, for determination of a stable isotope composition of hydrogen and oxygen, were stored in 50-mL HDPE bottles with double caps.

For chemical analysis, water samples were collected in dark glass bottles, triple rinsed with sample water prior to filling. For trace metal analysis, samples were stored in 50-mL HDPE bottles and acidified with suprapure HNO₃ to pH < 2.

2.4. Methods

2.4.1. Analyses

pH, electrical conductivity (EC), and water and air temperature were measured in situ in all spring and surface water samples. All chemical analyses were conducted in an accredited testing laboratory (under ISO 17025, International standard for testing and calibration laboratories) in the Company “Water supply and sewerage” in Podgorica. The coefficient of variation of hardness (CV), expressed in percent (%), is typically used to classify flow regimes in karst terrain [38]. Three flow types occur in a karst system: conduit, fissure, and diffuse flow (an equivalent to base flow) [39]. The coefficient of variation of hardness (CV) can be calculated by the following equation [40]:

$$\mathrm{CV}=100\frac{\mathsf{\sigma }\left(2.5{\mathrm{Ca}}^{2+}+4.1{\mathrm{Mg}}^{2+}\right)}{\overline{\mathrm{x}}\left(2.5{\mathrm{Ca}}^{2+}+4.1{\mathrm{Mg}}^{2+}\right)}$$

(1)

where Ca²⁺ and Mg²⁺ concentrations are in milligrams per liter (mg/L), the quantity in parentheses is the total hardness expressed as milligrams per liter CaCO₃, $\overline{\mathrm{x}}$ is the mean, and σ is the standard deviation of hardness.

The hardness of water is the concentration of ions in the water that will react with a sodium soap to precipitate an insoluble residue [41].

Stable isotope analyses of water were determined at the Department of Environmental Sciences of the Jožef Stefan Institute in Ljubljana (Slovenia) using the Finnigan MAT Delta plus isotope ratio mass spectrometer. Results are reported in the conventional delta (δ) notation (δ²H and δ¹⁸O), i.e., the relative deviation of the heavy-to-light isotope ratio of the sample from that of the standard (VSMOW) expressed in per mil (‰). In-house working standards calibrated vs. VSMOW2 and SLAP2 international reference materials were used to calibrate the measurements. The accuracy was checked using the USGS45- and USGS47-certified reference materials as controls, randomly distributed in each batch. The measurement uncertainty (determined as long-term deviation of control materials from their respective certified δ values) was 0.05‰ for δ¹⁸O and 0.7‰ for δ²H.

2.4.2. Calculation of Mean Transit Time (MTT) and Young Water Fraction (F_yw)

The assessment of the groundwater mean transit time (MTT) in the aquifer was based on the sinusoidal fluctuation of monthly stable isotope compositions (δ¹⁸O, δ²H) of precipitation and spring water. From the relationship between the amplitudes of the modeled sine wave, the MTT and young water fraction (F_yw) were calculated. The sine wave model has previously been used to fit seasonal variations of δ¹⁸O and δ²H in precipitation and spring water [42,43], defined as

δ = δ_mean + A(cos(c × t − θ)),

(2)

where δ is the modeled δ¹⁸O or δ²H, δ_mean is the mean annual measured δ¹⁸O or δ²H, A is the modeled annual amplitude of δ¹⁸O or δ²H, c is the radial frequency of annual fluctuations (0.017214 rad/day; [42]), t is the time in days after the start of the sampling period, and θ is the phase lag or time of the annual peak of δ¹⁸O or δ²H in radians [42,43].

For the purposes of estimating the MTT, the exponential model was applied [42,43], according to which

$$\mathrm{MTT}={\mathrm{c}}^{-1}{\left({\left(\frac{{\mathrm{A}}_{\mathrm{s}}}{{\mathrm{A}}_{\mathrm{P}}}\right)}^{-2}-1\right)}^{0.5}$$

(3)

where A_s and A_p are the modeled annual amplitudes of δ¹⁸O or δ²H of the spring water and of precipitation, respectively [43].

The young water fraction (F_yw) was estimated from the amplitude ratio, as suggested by Kirchner [44]:

$${\mathrm{F}}_{\mathrm{yw}}=\frac{{\mathrm{A}}_{\mathrm{s}}}{{\mathrm{A}}_{\mathrm{p}}},$$

(4)

3. Results and Discussion

3.1. Physicochemical Parameters

All results of the major physicochemical properties of the spring water (S1, S2, S3, and S4) and the Zeta River (RU and RD), during the research period, are summarized in Table S1 (Supplementary Materials), while the average, minimum, maximum, and median values are presented in Table S2 (Supplementary Materials). In general, all four sampled springs showed very similar temporal profiles of temperature, pH, EC, and HCO₃⁻ contents; also, the river water showed no significant difference regarding the measured parameters between the upstream and the downstream sampling sites (Figure 4).

The water temperature at all four springs was relatively uniform, with a mean value of 11.3 °C and annual variation between 6 and 13.3 °C, which is much lower than the local mean annual air temperature (around 16 °C). These results indicate a distant water infiltration area with lower mean annual air temperature and higher altitude [45]. The temperature of the river and the springs was very similar in the colder part of the year (November–May), with an average value of 10.8 °C. In the warmer period (June–October), the temperature of the Zeta River reached up to 21.3 °C (mean temperature of 17.8 °C), which is much higher than in the springs (mean temperature of 12.4 °C) (Figure 4a).

The pH value of the springs and the river showed little seasonal variation, ranging from 7.35 to 7.96 for the springs and 7.57 to 8.05 for the river. A significant difference in pH between the spring and river water (ranging from 0.2 to 0.4 units) was observed in the summer months when the river water temperature increased (June–September; Figure 4b). The increased pH value of the Zeta River in the warm period of the year may be due to a low water level, CO₂ degassing (see the discussion on pCO₂ below), and higher photosynthetic activity [46].

The temporal EC profiles (Figure 4c) of the spring and river water showed more variability than the other parameters. The variations of EC can help in estimating the relative residence time of the circulating water [47], and therefore, the response of the springs to rainstorms. During the dry season, similar trends of the EC in the river and the springs existed, while in the autumn and winter period, larger differences were recorded. In the periods of heavy rainfall, which caused turbidity of the spring water, the EC increased by approximately 20% (November 2017 and February 2018). This can be explained by the increased discharge of carbonate-rich groundwater masses from the aquifer. At the same time, the heavy rainfalls that caused turbidity of the Zeta River decreased its EC because of dilution with precipitation and surface runoff. In the rest of the winter period, the EC of the spring water was lower than that of the surface water (February, March, and December 2017, Figure 4c). The EC of the surface and spring water was evidently different during almost the entire research period, which indicates that the Mareza springs receive very little or no surface water from the Zeta River.

The coefficient of variation of hardness (CV) values for the Mareza spring water was 10% (calculated by using Equation (1)). According to the authors of [48], CV conduit flow-type springs varied by 10%–25%, while diffusive flow springs had a relative constant hardness and CV values below 5% (<5%). A conduit flow system is characterized by quick recharge and sensitive reaction after intense precipitation events, while a diffuse flow karst system is characterized by small variation of physiochemical parameters. The CV value of the Mareza spring (10%) is at the low end for a conduit flow system, indicating that the conduit system is influenced strongly by fractures that feed it or by diffuse recharge [40].

The content of the dominant anion (HCO₃⁻) in all samples of the Zeta River varied within a narrow range of 188–207 milligrams per liter, and in the spring water from 165 to 219 milligrams per liter (Table S2 and Figure 4d). A good positive correlation was found between bicarbonate and calcium ions of the spring water and the Zeta River water (r = 0.81) (Figure 5c), while no correlation was found between Mg²⁺ and bicarbonate. This indicates that the varying Ca²⁺ mainly originates from the dissolution of calcite.

In the study period, the Ca²⁺ between the sampling points along the Zeta River (RU and RD) differed by 0.1–11.79 milligrams per liter, with the maximum difference recorded in February 2018 (Figure 5a). The concentrations of Ca²⁺ and HCO₃⁻ fluctuated seasonally—in particular, in the period of heavy rainfall. Also, Ca²⁺ fluctuations were much greater than Mg²⁺ fluctuations. The concentration of potassium (K), in the greatest number of samples, was below the limit of detection (LoD; <0.5 milligrams per liter).

A useful parameter for better understanding the nature of the carbonate aquifer is the water–rock interaction through which the water circulated is a molar ratio of Mg/Ca (Figure 6), which depends on the proportion of calcite and dolomite present in the aquifer rock [39]. According to some authors, this parameter can be used as a qualitative indicator of the residence time of the water in the aquifer [49,50,51,52].

According to [39], the Mg/Ca molar ratio is 1 in water that dissolves pure dolomite, and 0 for water that dissolves pure calcite. The dissolution of calcite and dolomite in a 1:1 ratio results in an Mg/Ca ratio of water of 0.33. Here, all of the collected springs and surface water samples were below the line of 0.33 (Figure 6b), which suggests that limestone dissolution prevails over dolomite dissolution.

In spring samples collected in the period of heavy precipitation (14.11., 30.11., and 3.02.), we recorded a decrease of the Mg/Ca ratio (Figure 6a). At the same time, the spring water showed turbidity as a response of the karst spring to precipitation after an antecedent dry period. That turbidity events can explain the gradual increase in water level in the deep subterranean channels when, after the first heavy rain events, the hydraulic head within the karst system was high enough to mobilize the deep and long residence time water toward the surface as spring water. Considering the geological composition, the water turbidity was caused by the outwash of clay which occurs in cracks or pockets in karstified carbonate rocks [35].

The Mg/Ca time series of the Zeta River water was parallel to that of the Mareza springs, but consistently slightly lower (Figure 6a). The differences in the Mg/Ca ratios of the river water between the two sampling sites could be attributed to the discharge of small tributaries and the diffuse recharge downstream of the RU site. Most of the tributaries only flow seasonally, in the autumn and winter period, which is what caused the higher variability of the Mg/Ca ratio at the downstream site (RD) in the wetter period of the year. The fact that the Zeta River discharge at RD is almost double compared with that of the RU site supports this interpretation.

The calcite saturation index (SIc) and pCO₂ were determined using the computer program PHREEQC (version 3), developed by the United States Geological Survey (USGS) [53] which is part of the United States Government. The saturation index (SI) is defined as the log of the ratio between the ion activity product (IAP) of calcite in the water sample and that at equilibrium [54], and indicates the ability of water to precipitate (for supersaturated water with SIc > 0) or dissolve (for undersaturated water with SIc < 0) calcite [55].

During the research period, the SIc of the spring water shows variation from −0.23 to 0.27. In the period of a high amount of precipitation (November and December 2017, and February 2018) and a higher discharge of the springs, the water was undersaturated with respect to calcite (Sic < 0). The exception was the sample collected during the first heavy rain event (14 November 2017), which was supersaturated (SIc = 0.15). Then, probably, infiltrated rainwater mobilized the deeper groundwaters to discharge at springs. This water was warmer and its EC was higher (Figure 4). Small positive values of SIc (0.04–0.07), which indicate near-equilibrium conditions [56], were recorded in February, March, and April 2017 and January 2018. In the rest of the research period, the water was oversaturated with respect to calcite (SI = 0.12–0.27), with the highest SI value in June 2017.

The partial pressure of carbon dioxide in spring water fluctuated from 10^−2.67 to 10^−2.17 atm (atmosphere), (average, 10^−2.5; Figure 7a), while in the Zeta River water, it ranged from 10^−2.8 to 10^−2.4 atm (average, 10^−2.6; Figure 7b). All of these values are higher than the atmospheric partial pressure of CO₂ (10^−3.4 at 400 ppmv CO₂) and indicate degassing of CO₂ from both spring and river water throughout the year [54]. The estimated values for the Mareza springs were similar to or slightly below those considered normal in the groundwater of carbonate aquifers: 10^−2.5 atm [57]. The exceptions were the samples collected in November (30.11.), with the highest (maximum) value (10^−2.17 atm).

The sulphate concentration in the river water was low and showed little seasonal variation (2.74–5.19 mg/L), in contrast to the spring water, where sulphate content increased during the warmer and drier part of the year (Table S2 and Figure 8a). The molar ratio of SO₄/Cl vs. time (Figure 8b) showed a trend similar to sulphate concentration—the values were higher in spring water in the period with low precipitation and higher temperature, which supports the assumption of a geogenic rather than a meteoric source of sulphate.

In conclusion, the chemical composition of spring and surface water does not change much during the year, but nevertheless, some parameters show conspicuous patterns (e.g., electrical conductivity, calcite saturation indices). The most significant differences were noticed in the period of heavy rainfall, reflected in all major chemical parameters. Furthermore, all water samples showed the same water facies (calcium–hydrogen carbonate), low alkaline water, classified as low mineralization waters with low and insignificantly changing concentrations of chloride, sodium, and potassium during the research period. The biggest temporal variation was observed in the EC values closely related to the concentration of major ions (Ca²⁺ and HCO₃⁻), recorded in periods of heavy rainfall. During the research period, all of the springs (S1, S2, S3, and S4) showed the same chemical characteristics, while the surface water samples (RU, RD) showed increasing major ion contents downstream. According to chemical parameters, the main recharge area of the Mareza springs is located at a higher altitude (considering the constantly low temperature of the springs), with only a minor or negligible contribution from the surface water (Zeta River).

3.2. Isotope Composition of Oxygen and Hydrogen

For better understanding of the functioning of karst aquifers, the stable isotopes of water, oxygen, and hydrogen are most commonly used as tracers of water source, flow, and mixing, while the hydrochemical parameters provide information on the interactions of water with bedrock (e.g., dissolution and precipitation of minerals during infiltration and groundwater flow, and dilution with surface discharge). This combination of the stable isotopes of hydrogen and oxygen is also useful tools for determining residence time [42,58].

3.2.1. Isotopic Composition of Precipitation

At stations P1 and P2, 11 and te10n monthly precipitation samples were collected, respectively. Because of the dry weather, no samples were obtained at any of the sampling sites in October 2017, while in June 2017, enough samples could be retrieved only at station P2. No snowfall was recorded in the sampling area during the observation period, so all samples represented rainwater. The results of the δ²H and δ¹⁸O analyses are shown in Table S3.

Because of the very strong correlation between δ¹⁸O and δ²H values (r² = 0.93), further discussion on time series of stable isotopes will be based only on δ¹⁸O values.

The correlation between δ²H and δ¹⁸O in atmospheric precipitation, commonly called Global Meteoric Water Line (GMWL), is defined as: δ²H = 8δ¹⁸O + 10‰ [59], and the deuterium excess (d excess) parameter is defined as: d = δ²H − 8δ¹⁸O [33,60,61]. A deuterium excess value of 10‰ in the GMWL is based on global atmospheric water vapor that forms at a relative humidity of approximately 85%, producing a precipitation line that shifts from the seawater line by 10‰ [62]. Depending on the region, differences in the amount of precipitation, temperature variations, distinct air mass sources, evaporation, and fractionation processes occurring below the cloud base are characteristic at a local scale, which causes the relationship between the stable isotopes of water δ²H and δ¹⁸O to deviate from that of the GMWL [18]. The Local Meteoric Water Lines (LMWL) for particular regions can deviate from the GMWL, both regarding the slope and the d excess, depending upon the isotopic characteristics of precipitation at the local scale [61,62,63].

The calculated annual weighted average δ¹⁸O value of precipitation for the investigated period at station P2 was −5.9‰. The δ¹⁸O and δ²H values of precipitation show seasonal variations typical for continental stations in the Northern hemisphere, with more negative values in winter and considerably fewer negative values in summer [64] (−9.09‰ in February 2018 and −2.13‰ in July 2017; Figure 9). The isotopic compositions of O and H in precipitation in the study area are presented in Table S3 (Supplementary Materials). The most negative values at the station P2 were measured when the temperature was the lowest (average value of 7.4 °C) and the amount of precipitation was very high (total month amount 284 L/m²). Since these were the first analyses of the isotopic composition of precipitation in Montenegro, comparison with previous years was not possible.

Seasonal changes in temperature and amount of precipitation affect the isotopic composition of precipitation, which is reflected in the time series of δ¹⁸O in the precipitation collected at station P2 (Figure 9).

An observation period of one year is too short to estimate the Local Meteoric Water Line with reasonable uncertainty. For the purpose of this study, we could only estimate the short-term local meteoric trend line (δ²H = 6.94 δ¹⁸O + 5.87; Figure 9c), which may not be representative of long-term precipitation.

The slope and intercept of the trend line were lower than those of the GMWL, which suggests evaporation of the falling rain [65] as a result of the high temperature and the low humidity in the research area, in particular in the summer months.

The deuterium excess is commonly used as an indicator of water vapor source region [66,67], since it is sensitive to the meteorological conditions at the point where the vapor is originally evaporated from the surface, including sea surface temperature and relative humidity [68,69,70,71]. Monthly variations of d in the analyzed period are shown in Table S3 and Figure 9b.

Higher d excess in the autumn–winter period compared to spring–summer is typical of the Northern hemisphere [64]. In the Mediterranean area, precipitation has higher d excess than Atlantic air masses, generally increasing in W–E direction from 14‰ to 22‰ [72]. In the period from May to November 2017, at both measuring stations, the d values were lower than or close to the 10‰ that is typical for Atlantic air masses [73,74], and may be indicative of secondary evaporation processes (e.g., evaporation of falling raindrops) in a warm and dry atmosphere. Đorđević et al. [75] calculated the 96-hour backward trajectories ending at Herceg Novi, a coastal town 60 km West of Podgorica, and estimated that at more than 36% of days with 0.5 mm or more precipitation, the air masses derive from the western direction, originating in the northern Atlantic. This is in line with the findings of Schicker et al. [76], who estimated that in the Eastern Adriatic area, only a minor fraction of precipitation derives from the evaporation of local (Mediterranean) moisture. Therefore, the d excess values of precipitation in the study area typical for Atlantic air masses are not unexpected. In the winter period (November–May), these values range from 12‰ to 18‰, which are attributed to the precipitation originating from the Mediterranean Sea [72,74], with the exception of the results of March and July 2017, with a drastic deviation between two sampling locations (Table S3). This means that the recorded seasonal changes can be explained by the fact that the research area was under the influence of Mediterranean air masses during the winter period, while in the summer period, the research area was under the influence of air masses originating from the Atlantic Ocean.

The nearest station contributing data for long-term International Atomic Energy Agency (IAEA) Global Network of Isotopes in Precipitation (GNIP) is the Dubrovnik station, located 98 km WNW (West-Northwest) flight distance from the study area. The LMWL equation of Dubrovnik (δ²H (‰) = 6.46δ¹⁸O + 3.95) [77,78] for the period of September 2000–December 2003 is similar to the observed trend line for the study area (δ²H (‰) = 6.94δ¹⁸O + 5.9). However, because of the short observation period, the long-term d excess of precipitation in the investigated area may also considerably differ from the estimated value.

3.2.2. Isotope Composition of the Spring and River Water

All isotopic data of river and spring water samples are shown in Table S4, while the average, minimum, maximum, and median values are reported in Table S5 (Supplementary Materials).

During the study period, the δ¹⁸O values of the spring water ranged from −7.95‰ to −6.64‰ (Table S4 and Figure 10), and the δ²H values from −47.7‰ to −36‰ (Table S5); the ranges of values are much smaller than those in precipitation (see Table S3), suggesting good mixing of infiltrated water (precipitation) in the aquifer. All δ¹⁸O and δ²H values of the spring water samples plot between the GMWL and the Eastern Mediterranean MWL (EMMWL; [79]). Significantly lower δ values of the spring water compared to the local precipitation (Figure 9) point toward recharging with precipitation with a lower δ value than in the sampled area from a higher altitude. Assuming the altitude effect for the Southern Adriatic being about −0.26‰ per 100 m [77,80], the average altitude of the recharge area should be about 620 m asl. However, for a more reliable determination of the mean catchment altitude, it would be necessary to follow the long-term data of stable isotopes in precipitation and springs, including sampling stations at higher altitudes.

The most negative δ¹⁸O values of the spring water (down to −7.93‰) were recorded in the period between August and October 2017 (Figure 10), while the highest δ¹⁸O values (up to −6.7‰) were registered during the period with more precipitation (November–March). Meanwhile, the precipitation showed the most negative values in February 2018, while in the rest of the observation period, the δ values were much higher and showed no regular pattern. This indicates that the groundwater is homogenized in the aquifer, but no systematic shift between the isotope signal of spring water and precipitation in the study period could be detected. All spring water samples plot between the local precipitation trend line for the study period and the EEMWL, which leads to the conclusion that the origin of the spring water must be searched between the local precipitation and an area with precipitation closer to the EMMWL.

Figure 10 also shows the variation of the δ¹⁸O values of the river water between RU and RD (Table S4). The δ¹⁸O values measured at the downstream station (RD) are consistently higher than at the upstream station (RU), which is consistent with evaporation. This assumption is also supported by differences in the slope and intercept at both sampling points (δ²H vs. δ¹⁸O trend lines for RU: δ²H = 6.51δ¹⁸O + 3.48 and for RD: δ²H = 8.91δ¹⁸O + 22.1). The influence of several small ephemeral tributaries discharging to the Zeta River in the wetter part of the year between the RU and RD sites on the δ¹⁸O values of the river water cannot be ruled out; however, no isotope data for these streams are available.

3.2.3. Estimated Mean Transit Time (MTT) and Young Water Fraction (F_yw)

The data sets of the monthly δ¹⁸O and δ²H values in the precipitation and the spring water were used for the estimation of the MTT and F_yw of the Mareza aquifer.

The MTT was calculated using Equation (3), and the F_yw using Equation (4). The sine wave curves were previously modeled to fit the seasonal variations of stable isotope composition in precipitation and spring water (Figure 11 and Figure 12). Two separate analyses were performed, one with the δ¹⁸O data set and the other with the δ²H data set. The precipitation curves (Figure 11 and Figure 12) are based on the average values between the two precipitation sampling points (P1 and P2). Mareza spring S3 was selected as the most representative because it has the longest series (there are no significant differences in the isotopic composition between the four springs in Mareza). The radial frequency of annual fluctuations (c) is 0.017214 rad/day [43].

According to the estimation based on the δ¹⁸O data set, the MTT amounts 129 days, and the fraction of young water F_yw is 40.9%. The estimation based on δ²H values gave the following results: MTT = 92 days; F_yw = 53.3%. Thus, with two separate data sets, similar results were obtained. Since this is a first assessment from a rather short observation period, the results should be interpreted with caution.

4. Conclusions

The hydrochemical composition of the spring and river water showed seasonal variations, mostly depending on the amount of precipitation. Variable discharge influences the basic hydrochemical parameters, such as electroconductivity and the concentration of major ions (Ca²⁺ and HCO₃⁻). The springs generally showed little difference in hydrochemical composition, while in the river, the concentrations of the major solutes analyzed increased downstream between the two sampling sites. The recharge of the Mareza springs by the river water is possible, although the magnitude of such a recharge would be small. The isotopic composition of the spring and river water indicates that the main recharge area is in a mountainous area with conspicuously higher altitude than the investigated area.

The significant changes in hydrochemistry of spring water (EC, HCO₃⁻, Ca²⁺, and the saturation indices of calcite and pCO₂) recorded during the heavy rainfall lagged behind the change in discharge, suggesting that recharge water must first flush the system of stored water before arriving at the spring itself.

The isotopic composition of precipitation during the investigated period (2017–2018) in Podgorica and Danilovgrad was similar to that in Dubrovnik (Croatia)—GNIP station (at a distance of 98 km WNW).

The trend line for the δ²H vs. δ¹⁸O for precipitation in the investigated area was determined, although—due to a short observation period—it may have large attached uncertainties and therefore may differ from the long-term Local Meteoric Water Line. The seasonal variability of the δ¹⁸O values of precipitation and the d excess values can be attributed to the prevailing influence of Mediterranean air masses during the winter, while in the summer period, the investigated area was under the influence of air masses originating from the Atlantic Ocean.

The time series δ¹⁸O and δ²H values for the spring water (S1, S2, S3, and S4) showed the same behavior and the same origin of water. For precise determination of the mean altitude of the recharge area based on stable isotopes, a longer observation, and additional datasets on precipitation at high altitude meteorological stations would be necessary. The origin of the spring water must be searched between the local precipitation and an area with precipitation closer to the Eastern Mediterranean. The Zeta River, too, is recharged from water originating from both Mediterranean and Atlantic air masses, as indicated by the seasonal variability of the meteorological parameters. In the summer months, the evaporation of river water also influences the isotopic composition of the surface water.

According to the analysis of the δ¹⁸O and δ²H seasonal variations in the precipitation and the spring water, the groundwater mean transit time (MTT) of the Mareza aquifer is in the range of 92–129 days. Also, the young water fraction (F_yw) was estimated to be between 40.9% to 53.3%. The obtained values are typical of springs in highly karstified terrains [31]. Although these results are significant indicators of the degree of karstification and permeability of the Mareza aquifer, for now, they represent only the first assessment.

The results provide a good basis for future investigations and water management in the area, as well as in similar hydrogeological systems.