Multi-method geochemical characterization of groundwater from a hypogene karst system

An approach, combining several geochemical methods, was used to determine the groundwater properties and components of a hypogene karst system, where sampling is restricted only to the spring sites, and with a limited number of available sampling locations. Radiogenic isotopes (3H, 14C) were used to constrain the groundwater mean residence time and separate different groundwater components. Noble gases, stable isotopes of water (δ18O, δ2H), dissolved inorganic carbon (δ13C) and dissolved sulfate (δ34S, δ18O), and major ion and trace element composition were used to identify the source of water, its chemical evolution and water–rock interactions, as well as to identify the contribution and composition of endogenic gases. This approach was applied to three low-temperature thermal springs located in Mariovo (North Macedonia) associated with fossil hypogene caves, previously identified by morphological and geochemical studies of caves and cave deposits. Based on the obtained results, the main studied springs represent an output part of a regional hypogene karst groundwater system with a deep-circulating (~1 km), old (~15 ka), thermal (≥60 °C) water, which mixes with young (<50 years), cold (<14 °C) and shallow epigene karst groundwater. The output parts are structurally controlled, at the interception of low topography and deep faults, along which the groundwater interacts with deep-seated gases, dominantly CO2 of metamorphic origin (δ13C of +4.5‰ VPDB), with some contribution of mantle helium. The thermal karst groundwater interacts at depth with volcanic rocks from the nearby Neogene-Quaternary volcanic complex, as well as with metamorphic basement rocks and granitoids.


Introduction
Hypogene karstification is a relatively recent scientific paradigm in the understanding of the formation of karst systems (e.g. Klimchouk et al. 2017 and references therein). It addresses karst systems that develop due to bedrock dissolution by groundwater that recharges the karst rock formation from below (Klimchouk 2007). In contrast to the much more widespread epigene karst, where the chemical capacity to dissolve the bedrock is dominantly due to dissolved CO 2 acquired at the surface, more diverse genetic mechanisms are present in hypogene karst systems due to a variety of geochemical processes (Ford and Williams 2007;Klimchouk 2007;Palmer 2007).
As hypogene karst systems generally become accessible only after being intercepted by surface erosional processes, the most common approach to studying them is to characterize the cave morphology and deposits in fossil cave systems (e.g. Audra et al. 2002;De Waele et al. 2016;Plan et al. 2012;Spötl et al. 2016;Temovski et al. 2013). A complementary approach is to study the geochemical properties of the groundwater associated with hypogene karst systems.
Various geochemical methods and approaches have been used in the study of groundwater systems (Clark 2015). Many of them have been also applied to groundwater in carbonate aquifers that might be associated with hypogene karst processes (i.e. aquifers with thermal waters, high gas concentrations etc.), in order to identify sources of fluids, flow-paths and geochemical processes. A combined use of hydrochemistry and stable isotopes, including sulfur and strontium isotopes, has been applied to identify water-rock interactions and separate different groundwater flows associated with the thermal springs in Derbyshire, England, UK (Gunn et al. 2006). Gary and Sharp (2006) identified interaction of the groundwater at Sistema Zacatón (Mexico) with Pleistocene volcanic rocks based on strontium isotopes, and related the speleogenesis of the deep phreatic caves to volcanic derived H 2 S and CO 2 . Hydrochemistry and isotopic composition of dissolved inorganic carbon has been used to identify the source and amount of geogenic CO 2 related to the carbonate aquifers in the Appenine Region in Italy (e.g. Chiodini et al. 2000). A similar approach, combined also with noble gases, was used to identify geogenic gases in groundwater related to the giant collapse dolines in the Konya Closed Basin, Turkey (Bayari et al. 2009b(Bayari et al. , 2017, and sources of fluids in the deep phreatic cave Hranicka Abyss in Czechia (Sracek et al. 2019). Wynn et al. (2010) used sulfur isotopes to identify the sources of dissolved sulfur species in the thermal waters associated with sulfuric acid speleogenesis at Cerna Valley, Romania. Different radionuclides have been used to identify groundwater flow-paths and their mean residence times. Erőss et al. (2012) applied uranium, radium and radon to identify mixing of fluids in the Buda thermal karst system (Hungary) and to estimate the temperature and chemical composition of the end-members. Radiocarbon and tritium have been used to estimate the mean residence time of groundwater such as in the Bükk thermal karst system, Hungary (Hertelendi et al. 1995), the coastal thermal springs of Apulia, Italy (Santaloia et al. 2016), and the Konya Closed Basin in Tukey (Bayari et al. 2009a).
The approach in many of these examples involves using a relatively high number of sampling locations and/or time series of certain geochemical parameters covering longer periods (e.g. Bayari et al. 2009a, b;Erőss et al. 2020;Galdenzi et al. 2008;Mádl-Szőnyi and Tóth 2015). Spring waters often provide the main insight into the geochemical properties of these hypogene karst groundwater systems (e.g. Gunn et al. 2006;Mádl-Szőnyi and Tóth 2015;Palmer et al. 2017;Santaloia et al. 2016), due to their ease of access or due to lack of available boreholes. However, these discharge areas, generally located at points of low topography and controlled by major structural elements (Klimchouk 2007), are zones of convergence of discharge from different groundwater flow systems (local, intermediate and regional), allowing for mixing of groundwater with different geochemical composition (e.g. Erőss et al. 2012;Minissale et al. 2002). In such cases, to study the geochemical properties of a distinctive groundwater flow system, information about the geochemical composition of the end-members and their contributing ratios is required, which is not always attainable.
Here, an approach is presented of a combined use of several geochemical methods in the study of groundwater associated with hypogene karst systems, where sampling of the groundwater systems is restricted to only a small number of springs. Such an approach allowed for identification of the different components of groundwater where the end-member geochemical compositions were unknown. The approach was applied to Mariovo area in North Macedonia, where hypogene speleogenesis has been previously demonstrated based on morphological and geochemical studies of caves and cave deposits, and dissolution of carbonate rock due to cooling of CO 2 -rich thermal groundwater identified as the main speleogenetic process, coupled locally with sulfuric acid speleogenesis and ghost-rock weathering (Temovski et al. 2013;Temovski 2016). The insight obtained from this groundwater geochemical study was then used to provide an improved conceptual model of the hypogene karst system.

Research area
The study sites are located in the southern parts of the Republic of North Macedonia, in Mariovo, a hilly to mountainous area deeply incised by the valleys of Crna Reka and its tributaries (Fig. 1). The area spreads on both sides of the contact between the Vardar Zone and the Pelagonian Massif. These major tectonic zones are overlain by structures formed within the extensional tectonics of the South Balkan Extensional System (Burchfiel et al. 2008;Dumurdžanov et al. 2005), mainly the Mariovo Basin, bounded to the south by the mountainous terrain of the Kožuf-Kozjak volcanic centers (also known as Voras in Greece).
Karst is developed here in marble of supposedly Precambrian and Cambrian age, as well as in Cretaceous and Triassic limestone and dolomite, and varieties of Pliocene-Quaternary travertine. The Precambrian and Cambrian age of the Pelagonian marble formations has been recently questioned, assigning them Triassic to Jurassic age based on preliminary local 87 Sr/ 86 Sr measurements and various studies of its continuation to the south in Greece (Most 2003). The main structures are mainly spread in the NNW-SSE direction, with the carbonate rocks part of a series of nappe structures, except for the travertine deposits which are found either within, or topping, the Pliocene-Quaternary sedimentary sequences, or on Pleistocene river terraces. This long NNW-SSE karst stripe is cut transversely by the gorge-like valleys of Crna Reka (at Podot) and its tributary Buturica (at Melnica), where the three low-temperature (warm) thermal karst springs, that are the main study sites, are found ( Fig. 1).
Hypogene karstification in the area has been characterized on the basis of cave morphology and cave deposits (Temovski 2016). The main speleogenetic mechanism is hydrothermal speleogenesis, i.e. dissolution of carbonate rock due to cooling of CO 2 -rich thermal water, with increased geothermal gradient attributed to the Neogene-Quaternary Kožuf-Kozjak volcanism. At places, due to local geological or lithological control, this is coupled with other hypogenic processes such as sulfuric acid speleogenesis (Temovski et al. 2013 or ghost-rock weathering, i.e. preferential dissolution of calcite due to cooling of thermal waters, leaving in-situ dolomitic sand residue (Temovski 2016(Temovski , 2017. Melnica spring is located at the riverbed in the small gorge of Buturica Valley, with groundwater discharging from calcite marble along a lens of schists (Fig. 2). It is located 100 m below Provalata Cave, a fossil hypogene cave where both CO 2 -based hydrothermal and sulfuric acid speleogenesis were identified (Temovski et al. 2013. In Crna Reka, at Podot locality ( Figs. 1 and 3), the active thermal cave Karši Podot is developed in dolomite marble on the right valley side, and on a 3-m high river terrace on the opposite side of the river, the lukewarm Gugjakovo springs discharge from travertine deposits.
Preliminary findings from a previous sampling campaign in 2016  indicated that the water at these springs is a mixture of cold and thermal groundwater, with possible interaction with the Kožuf-Kozjak volcanics.
In addition to these, three more springs were also included: Manastir, Toplek 2 and Toplek 4. Although they do not belong to the same groundwater systems, they were added as comparative sites, representing either a young fresh groundwater system (e.g. Manastir and Toplek 4) or a low- Fig. 1 Location and geological setting of the study area. The approximate boundary between the Pelagonian Massif (PM) and the Vardar Zone (VZ) is also indicated. The potential recharge areas for different springs are also indicated: Gugjakovo Springs and Karši Podotcyan outline; Melnica Springcyan crosshatch; Manastirgrey crosshatch; Toplek 2 and Toplek 4black crosshatch temperature thermal groundwater system with similar geochemical composition (e.g. Toplek 2). Manastir is a cold spring, discharging from a large plateau of travertine deposits topping the Mariovo Pliocene-Pleistocene sedimentary sequence. Toplek 2 and Toplek 4 are located further east, on the opposite side of Kozjak Mt, in the foothills of Kožuf Mt (Fig. 1). Toplek 2 has similar temperature to Melnica and Karši Podot, while Toplek 4 is cold and intermittent. Toplek 2 and 4 are located very close to each other and discharge from dolomite rocks, near the Allchar Carlin-type ore complex (Palinkaš et al. 2018). Previous geochemical studies of Toplek 2 (Boev and Lepitkova 2003;Boev and Jančev 2014) show high concentrations of elements characteristic of the ore deposits (e.g. As, Sb, Tl).

Methodological approach
As the sampling network for the groundwater systems was restricted to only their springs (no available boreholes), and located in a remote area, with a limited number of sampling locations, the approach relied on the combined use of several geochemical methods.
Radiocarbon ( 14 C) and tritium ( 3 H) combined with noble gases ( 3 H-3 He method) were used to identify the mean residence time of the groundwater and to separate young and old components. Water stable isotopes and noble gases were used to identify the source of water and conditions at groundwater recharge. Major ion and trace element composition of the water, as well isotopic composition of the dissolved inorganic carbon (DIC) and dissolved sulfate, were used to interpret the chemical evolution of the groundwater and identify possible interaction with nonkarstic rocks. The concentration and carbon isotope composition of the DIC were used to identify different carbon sources and model the carbon isotope composition of endogenic CO 2 . The obtained fractions for different carbon sources (soil and endogenic CO 2 and carbonate rock) were then used to correct the radiocarbon age of the old groundwater by constraining the contribution of 14 C-free carbon. Noble gas concentrations and isotopic composition were also used to identify endogenic noble gases. The applied methodological approach is conceptually illustrated in Fig. 4, showing the use of characteristic geochemical parameters to identify different properties of the groundwater system. Further explanation is also given on the approaches and calculation steps used to identify CO 2 sources (section 'Sources of CO 2 and deconvolution of carbon components in DIC'), noble gas recharge temperatures (section 'Noble gas recharge temperatures'), different groundwater components, and their carbon isotopic composition and mean residence times (sections 'Calculating the age and fraction of the young groundwater component', 'Estimating the carbon isotope composition of the groundwater end-member', 'Correction for 14 Cfree dilution of DIC and calculating the age of the old groundwater').

Analytical methods
Field sampling was carried out in September 2018 at all locations, except Toplek 4, which was sampled in March 2019 when also Toplek 2 was resampled. Spring water temperature, electrical conductivity (EC) and pH were determined on-site using a Multi 340i multiparameter instrument. Total alkalinity was determined in the field by acidimetric titration, using HCl Samples were collected for laboratory analyses to determine water chemistry, stable isotopes, radiocarbon, tritium and dissolved noble gases. Water stable isotope composition, temperature, pH and EC of the three main springs were analyzed more frequently in addition to the main field campaign, when samples for full geochemical analysis were collected. The full list of those measurements is given in the electronic supplementary material (ESM).
Samples for determination of anions were filtered and stored without additional treatment, whereas samples for determination of cations and trace elements were acidified after filtration by addition of suprapure nitric acid. Major anion concentrations were measured on a Metrohm 850 Professional IC, and cations on an Agilent 4100 MP-AES. Major ion charge balance error was less than 5%. Trace element concentrations were measured using Agilent 8800 Triple Quadrupole ICP-MS. Stable isotope analyses of water and DIC were done on an automated GASBENCH II sample preparation device attached to a Thermo Finnigan Delta PLUS XP mass spectrometer (Vodila et al. 2011). Hydrogen and oxygen isotopes of the water are given as δ 2 H and δ 18 O values relative to Vienna Standard Mean Ocean Water (VSMOW), and DIC carbon isotopes are expressed as δ 13 C values relative to Vienna Pee-Dee Belemnite (VPDB). The precision of the measurements is better than ±0.15‰ for δ 18 O, ±2‰ for δ 2 H, and ± 0.08‰ for δ 13 C. Oxygen and sulfur stable isotope composition of the dissolved SO 4 2− were analyzed with a Thermo Finnigan Delta PLUS XP mass spectrometer attached to an elemental analyzer (Flash EA). Oxygen and sulfur isotopes of the dissolved SO 4 2− are given as δ 18 O values relative to VSMOW, and as δ 34 S values relative to Vienna Canyon Diablo Troilite (VCDT). An individual δ 18 O and δ 34 S measurement has a standard deviation (SD) of ±0.3‰ and ± 0.4‰, respectively.
Radiocarbon ( 14 C) analyses of the water samples were done by accelerator mass spectrometry (AMS), following standard procedures (Molnár et al. 2013a), on an EnvironMICADAS 14 C AMS facility (Molnár et al. 2013b). The DIC fraction of the water samples was extracted by acid (85% H 3 PO 4 ) reaction in a vacuum-tight reaction vessel. The produced CO 2 gas was cleaned and separated along an on-line gas handling system and graphitized by the sealed tube Znbased graphitization method (Rinyu et al. 2013). Radiocarbon results are expressed in absolute pMC units (percent Modern Carbon). The water samples for dissolved noble gases were Fig. 4 Conceptual representation of the applied methodological approach with indication of the analyzed geochemical parameters that were used to identify different properties of the karst groundwater system stored in copper tubes sealed by stainless-steel pinch-off clamps. Noble gas measurements were performed by a VG5400 (Fisons Instrument) and a Helix SFT (Thermo Scientific) noble gas mass spectrometers (Papp et al. 2012). Noble gas concentrations are given as cubic centimeters of gas at Standard Temperature and Pressure per gram (cm 3 STP/g). Tritium ( 3 H) concentration in water was determined using the 3 He ingrowth method (Palcsu et al. 2010). Tritium results are expressed as tritium units (TU), with precision of the measurements better than 0.1 TU.
All analytical measurements were done at the Isotope Climatology and Environmental Research Centre, Institute for Nuclear Research, Debrecen, Hungary.
The distribution of aqueous carbon species and saturation indices for aragonite, calcite and dolomite, as well as pCO 2 , were calculated using PHREEQC Version 3 software (Parkhurst and Appelo 2013). Statistical analyses of chemical parameters were done with IMB SPSS Statistics v.23 using normalized values.

Sources of CO 2 and deconvolution of carbon components in DIC
To identify the sources of CO 2 the applied approach first separates the amount of carbon contributed to the DIC by dissolution of carbonates (C rock ) from the external carbon (C ext ), i.e. carbon from CO 2 not accounted for by dissolution of carbonate rocks (Chiodini et al. 2000(Chiodini et al. , 2004Crossey et al. 2009). C rock can be calculated from the water chemistry data, as C rock = Ca 2+ + Mg 2+ − SO 4 2− , where the amount of Ca and Mg derived from dissolution of carbonate rocks is corrected for dissolution of sulfates (Chiodini et al. 2000). C ext is then calculated as C ext = DIC -C rock. The δ 13 C of the external carbon can be estimated from the following equation (Chiodini et al. 2000): (δ 13 C ext × C ext ) = (δ 13 C DIC × DIC) − (δ 13 C rock × C rock ). For δ 13 C rock the average values of measured samples of the carbonate rocks from the studied aquifers were used (Fig. 5).
The obtained δ 13 C ext ranged from values characteristic of soil-derived CO 2 , to values reflecting mixture of soil CO 2 and variable amounts of endogenic CO 2 . To estimate the δ 13 C of the endogenic CO 2 , the mixing of a soil end-member (with δ 13 C defined by samples with δ 13 C ext reflecting soil-derived CO 2 and a range of values for C ext ), with an endogenic endmember was modeled. The δ 13 C value of the endogenic endmember was constrained by the best fit mixing line between C ext and δ 13 C ext values of springs that are located next to each other, which, based on their geochemical composition and proximity, represent mixtures with variable contribution of the same end-members, and as such, they should fall on the same mixing line (section 'Stable isotope composition of DIC and the source of CO 2 '). The fractions of carbon in DIC derived from soil CO 2 (f soil ) and endogenic CO 2 (f endogenic ) at each spring can then be calculated from the binary mixing model using the equation: f soil = [(δ 13 C ext -δ 13 C endogenic )/ (δ 13 C soil -δ 13 C endogenic )]/ f ext , where δ 13 C soil is the soil endmember value, δ 13 C endo is the modeled endogenic endmember value, f ext is the fraction of external carbon in DIC calculated as f ext = C ext /DIC and f endogenic = 1 -f soil . The carbon in DIC derived from dissolution of carbonate rocks (f rock ) can be calculated as f rock = C rock /DIC, and the total carbon balance is expressed as DIC = f soil + f endogenic + f rock = 1.

Noble gas recharge temperatures
The noble gas recharge temperatures (NGTs) reflect the temperature during recharge at which the groundwater was equilibrated to the soil air noble gas concentrations (Aeschbach-Hertig and Solomon 2013). The NGTs were estimated by using all noble gases except He, which usually has higher concentration due to contribution from other sources (crustal, mantle). The noble gas concentration in air, besides the temperature, depends also on air pressure, which varies with elevation. An estimation on the recharge elevation was made by extracting the mean elevations from a digital elevation model (DEM) for the outcrops of the carbonate formations likely contributing to these aquifers ( Fig. 1). NGT values were calculated using the NobleBook excel worksheet (Aeschbach-Hertig et al. 2000. As the recharge elevations are an average value, the estimated NGTs are also reflecting an average recharge temperature of water infiltrated at supposedly different elevations. The uncertainty in the NGTs due to recharge elevation uncertainty (taken as 1SD from mean elevations) ranged from 0.4 to 1.6°C, which is 4-11% of the calculated NGT value.
Calculating the age and fraction of the young groundwater component The age of the young groundwater cannot be determined based only on 3 H data, due to possible 3 H contribution from the thermonuclear testing in the 1950-1960s (bomb peak 3 H). By combining 3 H and noble gas data, based on the accumulation of 3 He ( 3 He trit ) from the decay of 3 H (half-life of 12.32 years), an apparent ( 3 H-3 He) age (expressed in years) can be calculated using the decay formula t = −17.77 × ln[ 3 H/ ( 3 H + 3 He trit )], which represents the period since the 3 H-bearing water reached the groundwater (Schlosser et al. 1988). 3 H in the water recharging the aquifer before the thermonuclear testing, i.e. older than 1953 (assuming 4 TU for 3 H in precipitation), would have decayed by the sampling time (2018) to a value of <0.1 TU. When plotted on a graph showing the 3 H in precipitation decayed to the sampling time, samples with >0.1 TU that are below the precipitation curve represent mixture of old 3 H-free groundwater and young groundwater. By using the 3 H-3 He apparent ages to constrain the recharge period of the young component, the fraction of the young groundwater can be estimated from the measured 3 H concentration of the groundwater and the expected recharge 3 H concentration obtained from the precipitation 3 H dataset for the given recharge period, decayed to the sampling time. For the precipitation 3 H concentrations, as there are no historical data available for the area, the data from the Vienna GNIP (Global Network of Isotopes in Precipitation) station was used, both monthly and yearly datasets.

Estimating the carbon isotope composition of the groundwater end-members
If it is assumed that the young groundwater DIC evolved due to carbonate dissolution under a system closed to soil CO 2 , then the 14 C of the young component can be calculated from 14 C young = 0.5 × ( 14 C CO2 + 14 C rock ) (Han and Plummer 2016), i.e. half of the 14 C of the soil equilibrated CO 2 , as 14 C rock is 0 pMC. The 14 C of the soil equilibrated CO 2 can be estimated from the 14 C of the atmospheric CO 2 (Hua et al. 2013) for the given recharge period based on the 3 H-3 He age. By using springs that represent mixtures with variable contribution of the same end-members for the young and old groundwater (e.g. Gugjakovo and Karši Podot in this case), the δ 13 C of the young component can be calculated from the linear regression of their δ 13 C and 14 C values ( 14 C = a + b × δ 13 C) and the known 1 4 C value for the young component (i.e. δ 13 C young = ( 14 C young -a)/b; Fig. 6). The carbon isotope composition (δ 13 C and 14 C) of the old groundwater can then be calculated from the composition of one of the springs and the young groundwater by using a binary mixing model [C mix = f young × C young + (1 -f young ) × C old ] and the fraction of the young component (f young ) for selected spring estimated from 3 H (Fig. 6).
Correction for 14 C-free dilution of DIC and calculating the age of the old groundwater The radiocarbon content of the DIC can be used to estimate the mean residence time of the old groundwater component, assuming that the lowering of DIC 14 C content along the groundwater flow-path is only due to 14 C decay (half-life of 5,730 years). However, there can be a number of processes and mechanisms along the flow-path (e.g. carbonate rock dissolution, methanogenesis, endogenic CO 2 , mixing with young groundwater etc.) that can alter the 14 C content of DIC, masking the decay signal (Clark 2015;Han and Plummer 2016). In epigene karst groundwater systems, carbonate dissolution is the major modifying factor, diluting the soil 14 C signature by addition of 14 C-free (dead) carbon, as this is the main process that controls the formation of the system (Ford and Williams 2007). In hypogene karst systems, this can be further complicated by the presence of various endogenic gases, that can contribute dead carbon to the DIC pool, especially CO 2 of metamorphic or magmatic origin (e.g. Bayari et al. 2009a;Clark et al. 1989). To correct for the 14 C dilution, a number of methods and approaches have been developed (see Han and Plummer 2016 for a recent review). Most of them are based on the mass balance of DIC carbon species or carbon isotopes and determine a radiocarbon dilution factor q, that represents the fraction of the initial 14 C at recharge (e.g. Ingerson and Pearson 1964;Tamers 1975;Fontes and Garnier 1979;Mook 1980). Their applicability is constrained by the assumed chemical reactions and isotopic fractionation along the evolution of the DIC carbon composition (Han and Plummer 2016). Some of them are based on a statistical approach and rely on the use of a larger number of samples from the same groundwater system (e.g. Gonfiantini and Zuppi 2003), or they apply geochemical modeling for the evolution of 14 C along the flow-path of the groundwater system (e.g. Plummer et al. 1994;Bayari et al. 2009a). Mixing with younger groundwater can further complicate the 14 C signal, by contributing DIC with higher 14 C content due to the thermonuclear testing in the 1950-1960s; thus, in some studies (e.g. Gonfiantini and Zuppi 2003), groundwater samples that indicate mixing with younger (bomb-peak affected) groundwater, are excluded from the determination of the 14 C-based mean residence time. However, sometimes the mixed groundwater is the only available option, especially when sampling is limited to spring sites, and can still provide valuable information about the groundwater system.
For the selected study area, the DIC 14 C was altered by addition of dead-carbon from dissolution of carbonate rocks, and by endogenic CO 2 , but also by mixing with younger groundwater (section 'Radiocarbon and the mean residence time of the old groundwater'). To correct for this, the used approach was based on the previously deconvoluted components of carbon derived from dissolution of soil CO 2 , carbonate rock and endogenic CO 2 (section 'Sources of CO 2 and deconvolution of carbon components in DIC'). The soil carbon component of the old groundwater (f soil-old ) can be calculated from the young groundwater fraction (f young ), given that the endogenic component is provided solely by the old groundwater and the young groundwater has half of the carbon from the bedrock and half from the soil (closed system carbonate dissolution), i.e. f old = (f rock -0.5f young ) + (f soil -0.5f young ) + f endogenic . The soil fraction in the old groundwater is then calculated from f soil-old = (f soil -0.5f young )/f old , and represents a different derivation of the dilution factor q, that also considers contribution of dead carbon from endogenic CO 2 . The radiocarbon age of the old component (expressed in years) can then be calculated from the decay equation t = −8267 × ln[ 14 C old /(q × 14 C 0 ], where 14 C old is the old groundwater 14 C composition and 14 C 0 is the initial soil 14 C at recharge, taken as 100 pMC.
HCO 3 is the dominant DIC component, with higher CO 2 concentrations found at Gugjakovo, Karši Podot and Melnica. Calculated pCO 2 values at Manastir, Toplek 2 and Toplek 4 are lower than 10 -1.7 atm (Table 2) and within the range for soil derived CO 2 (10 -2.5 to 10 -1.5 atm; Clark 2015). pCO 2 at Melnica and Karši Podot (~10 -0.6 atm) is higher than expected Toplek 4 (0.36) reflect their dolomite aquifers lithology, with higher Mg/Ca ratios at Toplek 2 likely reflecting the longer water-rock interaction. The studied springs have higher concentrations of certain trace elements (TE), of which boron, strontium, arsenic, and lithium are the most abundant (Table 3; (Boev and Lepitkova 2003;Boev and Jančev 2014). Toplek 2 is located next to the Allchar As-Sb-Tl-rich ore deposit (Palinkaš et al. 2018), where the trace element composition was attributed to higher water-rock interaction with the volcanic rocks of the Kožuf-Kozjak volcanism. Lithium is increased at Melnica, Karši Podot and Gugjakovo (35-51 μg/L), with much lower concentrations (<3 μg/L) found at the other springs. Higher concentrations of boron and lithium found at Melnica, Karši Podot and Gugjakovo, but not at the Toplek springs, might be due to contact at depth with the granitoid plutons of the Pelagonian massif. Strontium is generally increased at all springs, except Toplek 4, with concentrations ranging from 82 to 348 μg/L, and Melnica having highest concentration.
Pearson's correlation coefficients of the physico-chemical parameters show several highly correlated parameters (Table 4; ESM), which reflect water-rock interactions with three main aquifer lithologies. HCO 3 , Ca, and Mg are reflecting the dominant dissolution of carbonate rocks. Na and K likely reflect water-rock interactions with the various schist formations or the Pelagonian metamorphic basement. Strong correlation of B, U, Cs, Rb and Sr points to contact with the granitoid bodies within the Pelagonian basement. The highly correlated As, Mo, Sb and Tl clearly reflect contact with the Kožuf-Kozjak volcanic complex or ore deposits related with earlier hydrothermal system connected to the volcanism. This is further illustrated by the hierarchical cluster analysis (HCA) which reflects the similarities between the studied groundwater (Fig. 8). The clustering of the chemical parameters reveals distinct groups, and in addition to the circulation through the carbonate formations, indicates also interaction with rocks from the Pelagonian basement, and rocks of the Neogene-Quaternary Kožuf-Kozjak volcanic complex.   Fig. 9). This can be a result of dominant recharge of the aquifers during the cooler period of the year or higher contribution from water infiltrated at higher elevation, or, for the older groundwater, it might indicate recharge under cooler climate conditions.

Stable isotope composition of the dissolved sulfate
The sulfur stable isotope composition of the dissolved SO 4 2− was analysed with attempts also to measure the oxygen stable isotopes. Melnica and Karši Podot had similar values for δ 34 S (−1.3 and −1.4‰), Toplek 2 had lowest values (−4.0‰), while positive δ 34 S values were found at Manastir (+2.1‰) and Gugjakovo (+5.1‰). δ 18 O values were positive and ranged from +0.2‰ at Toplek 2 up to +5.6‰ at Karši Podot (Table 5). The sulfate stable isotopes at Melnica, Karši Podot and Toplek 2 are within the range of values (−7.5 to +0.7‰ δ 34 S; −3.9 to +8.2‰ δ 18 O) found in the gypsum deposits of Provalata Cave , whose origin is attributed to oxidation of H 2 S in vadose settings (i.e. sulfuric acid speleogenesis by condensation corrosion). Their δ 34 S values are also within the range (−9.8 to 0.9‰) of the sulfide minerals found in Allchar ore deposit, where contribution of magmatic H 2 S was identified for the main mineralization stage (Palinkaš et al. 2018). The source of the dissolved sulfate at these springs, thus can be attributed to oxidation of deepseated sulfides (either H 2 S or sulfide minerals deposited by earlier hydrothermal systems). The values at Gugjakovo and Manastir indicate a different source, likely oxidation of sedimentary pyrite.

Stable isotope composition of DIC and the source of CO 2
The δ 13 C DIC showed clear differences between the springs (Table 6) The pCO 2 , pH and δ 13 C DIC at Manastir and Toplek 2 indicate DIC evolution by carbonate dissolution mostly under closed system conditions after some initial open system dissolution. The very low pCO 2 (10 -3.5 atm), high pH and low alkalinity at Toplek 4 indicate closed system carbonate dissolution. The high pCO 2, low pH and high δ 13 C at Melnica and Karši Podot indicate a possible contribution of endogenic CO 2 .
The obtained values for δ 13 C ext range from −20.8 to −2.3‰ (Table 6). Values reflecting soil δ 13 C CO2 are found at Manastir and Toplek 4 (−20.1 ± 0.7‰), reflecting the dominantly C 3 vegetation of the area (mostly covered by forests), and the others show variable contribution of endogenic CO 2 . To estimate the δ 13 C of the endogenic CO 2 , the mixing of a soil endmember with δ 13 C of −20.1‰ and a range of 0.5-5 mmol/L for C ext, with an endogenic end-member was modeled (Fig. 6). The C ext and δ 13 C ext values of Karši Podot and Gugjakovo were used to constrain the δ 13 C of the endogenic end-member, as their close proximity and similar geochemical composition indicate that they share the same end-members. The closest  --match was found for a mixing line defined by soil CO 2 with δ 13 C of −20.1‰ and C ext of 4.9 mmol/L to an endogenic endmember with δ 13 C value of +4.5‰ (Fig. 10). Such δ 13 C value indicates a metamorphic origin for the endogenic CO 2 component (Chiodini et al. 2000;Clark 2015). However, based on the noble gas data (section 'Noble gases, tritium and the mean residence time of the young groundwater'), there is some contribution of mantle helium, indicating also possible contribution of mantle CO 2 , thus the metamorphic CO 2 δ 13 C might be somewhat underestimated. Based on the three identified sources of carbon (carbonate rock, soil CO 2 and endogenic CO 2 ), the fractions of endogenic CO 2 range from 13% at Toplek 2 and 19% at Gugjakovo, to 51 and 54% at Karši Podot and Melnica, respectively.
Noble gases, tritium and the mean residence time of the young groundwater Noble gas concentrations are close to air-equilibrated values, except for He concentrations, which are somewhat higher at Melnica, and much higher at Gugjakovo ( Table 7). The R/R a ratios (where R and R a are the 3 He/ 4 He ratio of the sample and air, respectively), range between 0.76 and 1.16, with Toplek 4 reflecting air composition (R/R a ≈ 1), and the others having R/ R a values either lower or higher than 1.
The NGT values range from 3.8 to 22.7°C (Table 7). The estimated NGT for Manastir (11.8°C) is within the range of MAAT for the area (Fig. 11a). The NGTs for Gugjakovo, Toplek 2 and Toplek 4 are somewhat lower than expected for the local MAAT at the spring sites, which might be due to different reasons such as infiltration during older cooler periods, larger contribution of water infiltrated at higher elevation, or fast and/or higher contribution from water infiltrated during the winter period. The coldest estimated NGT at Toplek 4 (3.8°C), clearly reflects winter infiltration, also confirmed by the low spring water temperature, and the recent 3 H-3 He age. Their δ 18 O and NGT values follow the relationship of local mean monthly temperatures and monthly precipitation δ 18 O (as observed during 2018-2019; Fig. 11b), which might be a result of higher contribution of winter infiltration.  ), soil CO 2 (C soil ) and endogenic CO 2 (C endo ). Sample codes same as in Fig. 7 The estimated NGTs for Karši Podot and Melnica are much higher and clearly reflect the temperature of the groundwater at the spring (Fig. 11a). At Karši Podot, this is probably due to re-equilibration of the dissolved noble gases to the cave atmosphere behind the point where the water emerges at the cave wall. Considering the sponge-like morphology of the cave (Temovski 2016), there is likely a large contact of air with the water table along many small interconnected cavities. Similarly, at Melnica, a close match of the NGT to the spring water temperature might indicate an existence of a cave passage behind the spring where the water table is in contact with a large cave air mass. Toplek 4 (5.8 TU) had the highest 3 H values. Assuming no mantle helium, and after removal of air-derived and terrigenic (crustal) He, the calculated concentrations of 3 He formed by 3 H decay ( 3 He trit ), except for Melnica and Gugjakovo, range from 0 to 8.3 × 10 −14 cm 3 STP/g, yielding 3 H-3 He apparent ages of 0 ± 3 years for Toplek 4, 10 ± 3 years for Manastir, 31 ± 2 years for Karši Podot, and 49 ± 1 years for Toplek 2. The calculated 3 He trit concentrations were much higher at Melnica (68 × 10 −14 cm 3 STP/g) and especially at Gugjakovo (909 × 10 −14 cm 3 STP/g), indicating presence of mantle derived 3 He, thus no 3 H-3 He apparent ages were calculated.
The possible mantle helium contribution can be estimated from a three-component (atmospheric, crustal and mantle helium) mixing model, using R/R a values corrected for excess Table 7 Noble gas concentrations (cm 3 STP/g), helium isotope ratios and noble gas recharge temperatures (NGT). R c /R a is the entrapped air-and tritiogenic 3 He-corrected ratio air and tritiogenic 3 He (R c /R a ), and the 4 He/ 20 Ne ratio. For M e l n i c a a n d G u g j a k o v o , c o n s i d e r i n g t h e h i g h nonatmospheric He concentration, the 3 He trit contribution is small. Nevertheless, the 3 He concentration was corrected based on a 3 He trit value calculated using the maximum possible 3 H input value that can decay to the measured one in the period since the beginning of the thermonuclear testing (taken as 1953). Using 4 He/ 20 Ne ratio of 1000 and R/R a of 8 and 0.02 for the mantle and crustal end-member (Ozima and Podosek 2002), respectively, the calculated mantle He contribution at Melnica is 10% and at Gugjakovo 9%, with also small (2%) mantle contribution found at Karši Podot. Manastir and Toplek 2 have a mixture of crustal and atmospheric helium, while Toplek 4 has only atmospheric helium (Fig. 12).
The 3 H values of studied springs are lower than the decayed yearly 3 H values for Vienna precipitation (Fig. 13). As pointed out, this likely reflects the proportion of the young water fraction, although the actual magnitude could differ due to slightly different values for the local precipitation 3 H compared to the used dataset. The mixing of young and old water and their variable fractions are also indicated by other parameters such as the combined presence of 3 H with very low values for 14 C, and water temperature. However, it should be noted here that the 3 H-3 He age at Karši Podot might be an overestimation due to contribution of mantle 3 He; this is because 2% mantle helium was estimated based on the combined helium and neon isotope data. By using the 3 H-3 He apparent ages to constrain the recharge period of the young component, for Karši Podot, from the yearly precipitation 3 H dataset, 30-47% of young groundwater is found (22.4-128% if the monthly 3 H dataset is used), with the maximum value clearly an overestimation (Table 8). For Manastir, the estimated young groundwater contribution is 32-54% based on the yearly data, or 25-80% based on the monthly data. However, this might be an underestimation.
Although there are no local data available for the estimated recharge period, recent monitoring (since 2018) of 3 H in the local precipitation shows 3 H values as low as 4-5 TU for the winter period. If similar values are assumed for the estimated recharge period (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012), then the 3 H value at Manastir can be explained by 3 H decay from precipitation infiltrated dominantly during the winter period, thus the water represents only young groundwater. Toplek 2 is much lower than the precipitation curve, based on which the estimated young groundwater contribution is 13-19% (yearly data) or 10-35% (monthly data).

Radiocarbon and the mean residence time of the old groundwater
The measured DIC 14 C concentrations ( 14 C DIC ) range from 9.7 to 71.4 pMC, and also show three groups of values: Melnica (9.7 pMC) and Karši Podot (15.3 pMC) have the lowest, Toplek 2 (48.5 pMC) and Gugjakovo (48.9 pMC) intermediate, and Toplek 4 (68.4) and Manastir (71.4) the highest concentrations (Table 9; Fig. 14a). Based on their 3 H-3 He age and carbon isotope composition, Manastir and Toplek 4 are considered as modern waters, with their 14 C values diluted by carbonate dissolution under a system mostly closed to soil CO 2 . The waters at Melnica, Karši Podot and Gugjakovo represent a mixture of young and old groundwater, and from their geochemical composition it is evident that Gugjakovo has the highest and Melnica the lowest contribution of young groundwater.
As discussed before, Karši Podot and Gugjakovo likely represent mixtures with variable contributions of the same end-members for the young and old groundwater. Furthermore, from the 3 H-3 He age, the recharge period for the young component (1986)(1987)(1988)(1989) is known, for which the 14 C of the atmospheric CO 2 (Hua et al. 2013) was between Fig. 12 Helium components in the studied groundwater, corrected for excess air and tritiogenic 3 He: a three-component (atmospheric, mantle, crustal) helium mixing plot; b calculated percentages of mantle-, crustal-and air-derived helium. Sample codes (a) same as in Fig. 7 detailed in section 'Estimating the carbon isotope composition of the groundwater end-member'', the calculated carbon isotope composition of the young component is 58.8 pMC 14 C and − 8.9‰ δ 13 C (Fig. 14b). The old groundwater at Karši Podot and Gugjakovo (Podot locality), calculated from the carbon isotope composition of the young groundwater and Karši Podot and the estimated young fraction at Karši Podot, thus has 14 C of 2.7 pMC and δ 13 C of 0.07‰, reflecting the large contribution of endogenic CO 2 . From the binary mixing model and the estimated compositions of the end-members, 82% contribution of the young groundwater can be found at Gugjakovo.
To calculate a radiocarbon age for the old component at Podot (Gugjakovo and Karši Podot), first a correction has to be made for the dilution of the soil 14 C signature by addition of 14 C-free (dead) carbon from dissolution of carbonate rock and endogenic CO 2 . Using the deconvoluted components of carbon in the old groundwater (section 'Stable isotope composition of DIC and the source of CO 2 '; Table 6), and the approach detailed in section 'Correction for 14 C-free dilution of DIC and calculating the age of the old groundwater', a dilution factor of 0.17 can be calculated, from which a conventional 14 C age of 15.3 ka is obtained for the old groundwater at Karši Podot and Gugjakovo (Table 9).
For Melnica there is no separate estimate for the young fraction, so the same approach cannot be used to calculate the carbon isotope composition of the old component. However, the Melnica sample is very close to the δ 13 C-14 C mixing line of Podot locality (Fig. 14b), indicating that the old component at Melnica will not be significantly different. If no young groundwater contribution is assumed, the dilution factor q will be represented by the total of soil derived carbon (0.20), and from the measured 14 C of 9.7 pMC a 14 C age of 6.1 ka can be calculated. This should be considered as a minimum age, as young groundwater contribution is clearly indicated by the presence of 3 H. The 14 C age of 15.3 ka obtained for Podot locality thus should be considered as a maximum age for Melnica. Using the carbon composition of the young end-member at Podot locality, a maximum contribution of 12.8% young groundwater can be calculated for Melnica, as well as an old end-member composition of 2.5 pMC 14 C and 0.24‰ δ 13 C (Table 9).  Based on the 3 H data Toplek 2 also appears to be a mixture of old groundwater with 10-35% young groundwater. The δ 13 C of −9.3‰ at Toplek 4 can be considered as a representative value for the young component evolved under closed system carbonate dissolution. From this value and the fractions for the young component, a δ 13 C value of −7.4 ± 0.3‰ can be calculated for the old groundwater. The 14 C composition of the young component at Toplek 2, should be higher than the one of Toplek 4 as, based on the 3 H-3 He age, it must have been affected by bomb-derived 14 C. For the estimated recharge period (1969)(1970)(1971) the atmospheric 14 C was 151-156 pMC (mean 153.5 pMC), and with half of that (closed system dissolution) as the young component 14 C composition, the old groundwater 14 C at Toplek 2 should have 38.7 ± 6.8 pMC. The dilution factor q at Toplek, calculated from the carbon components is 0.46, using which a 14 C age of 1.5 ± 1.5 ka can be calculated (Table 9). The large age error at Toplek 2 is due to the large error in the estimate for the contribution of the young groundwater.

Geothermometry
An attempt was made to estimate the reservoir temperatures of the warmest springs using different geothermometers based on: SiO 2 (Fournier 1977), Na-K-Ca (Fournier and Truesdell 1973) with corrections for Mg (Fournier and Potter 1979) and CO 2 (Paceš 1975), as well as K-Mg (Giggenbach 1988) and δ 18 O SO4-H2O (McKenzie and Truesdell 1977). The estimated reservoir temperatures range from 60 to 175°C (Table 10). The most comparable temperature estimates are given by the SiO 2 , K-Mg and Na-K-Ca (pCO 2 corrected) geothermometers, with values of 60-65°C for Karši Podot, 60-71°C for Melnica and 64-88°C for Toplek 2. These estimates are considered here as minimum temperatures, as the waters were diluted by mixing with 10-20% cold water, although the cold water has likely much lower dissolved content. Palinkaš et al. (2018) estimated temperatures from~100°C to more than 200°C for the hydrothermal fluids that carried the mineralization of Au, As, Sb and Tl at Allchar. Considering this, 100°C can be considered as an upper limit for the reservoir temperature of these groundwater systems. The lower estimated reservoir temperatures than the earlier hydrothermal systems likely reflect the later stage of the thermal evolution of these systems. This is also supported by the estimated temperature of <50°C for the older phase (>1.5 Ma) of calcite deposition in Provalata Cave (located above Melnica Spring), that acted as a feeder to a former spring (Temovski et al. 2013;Temovski 2016).  Geothermal gradients are not well constrained, with reported values of more than 40°C/km for the Vardar Zone (Kotevski 1987). Using a geothermal gradient estimate of 40°C/km, the circulation depths are ranging from 0.9-1 km at Karši Podot, 0.9-1.2 km at Melnica and 1.1-1.7 km at Toplek 2. However, the geothermal gradient in this area is also likely underestimated (especially at Toplek locality), thus both the reservoir temperature and the circulation depth are considered as a loose approximation.

Discussion
Estimation of groundwater component fractions by combined use of 14 C, δ 13 C, 3 H and noble gases As indicated by several parameters in their geochemical composition (e.g. temperature (T), EC, TDS, 3 H, 14 C, δ 13 C), the studied spring waters at Melnica, Karši Podot and Gugjakovo represent mixtures of two end-members, a young fresh groundwater and an old thermal groundwater. Without clear end-member representatives at the studied locations, the approach to estimate their compositions is based on one main assumption: the selected springs (Karši Podot and Gugjakovo) share the same end-members, i.e. they are mixtures of different proportions from the same groundwater end-members. As indicated by their close proximity and especially their geochemical composition, this assumption seems valid for the two selected springs, with Karši Podot having higher contribution of the old groundwater (e.g. higher T, EC, TDS, δ 13 C and trace element concentrations, and lower 3 H and 14 C) and Gugjakovo having higher contribution of the young groundwater (e.g. lower T, EC, TDS, δ 13 C and trace element concentrations, and higher 3 H and 14 C). The fraction of the young component was estimated based on combined use of the spring water 3 H concentration and 3 H-3 He apparent age (young component recharge period), by comparison to the historical record of 3 H concentration in precipitation decayed to the sampling date. Assuming carbonate rock dissolution under a system closed to soil CO 2 , the carbon isotopic composition of the young component was estimated from atmospheric CO 2 ( 14 C) and the linear relationship of the measured carbon isotope composition at Gugjakovo and Karši Podot (δ 13 C). To test the reliability of this estimate, the δ 13 C of the soil equilibrated CO 2 from δ 13 C young = 0.5 × (δ 13 C CO2 + δ 13 C rock ), assuming closed system carbonate dissolution (Han and Plummer 2016). Using the δ 13 C value of the young component (−8.9‰) and a value of 2.8‰ for δ 13 C rock (section 'Sources of CO 2 and deconvolution of carbon components in DIC'), the calculated δ 13 C of the soil equilibrated CO 2 of the young groundwater is −20.7‰. This is close to the values obtained for the soil-derived external CO 2 at Manastir and Toplek 4, where closed system carbonate dissolution was found (section 'Stable isotope composition of DIC and the source of CO 2 ').

The source of the groundwater
Water stable isotopes, as well as NGTs, clearly show that the groundwater is of meteoric origin. The δ 18 O and NGT values reflect recharge under cooler temperatures than local MAAT, which, especially for the young shallow groundwater, probably reflects higher infiltration of water during the winter period. This is to be expected in temperate areas, where even though the precipitation amount might be lower in the winter period, higher evapotranspiration in the summer period prevents larger infiltration (e.g. Jasechko et al. 2014). However, at Karši Podot and Melnica, where the contribution of shallow young groundwater is small, even though the NGTs are not reflecting recharge values, lower δ 18 O and δ 2 H values might reflect lower MAAT values, which could be as a result of recharge of the old groundwater at cooler climate conditions, or considering the hydrogeological settings, recharge at higher elevations, or both.

Water-rock interactions
Although the trace element composition indicates interaction with volcanic rocks of Kožuf-Kozjak volcanic system, as well as metamorphic and magmatic rocks of the Pelagonian basement, the major ion composition clearly reflects the carbonate aquifer lithology. Mg/Ca molar ratios indicate dissolution of both calcite and dolomite minerals at Melnica, Karši Podot and Gugjakovo springs. The Mg/Ca ratio at Melnica, discharging from calcite marble, indicates also flow through the dolomite marble formation, which is quite plausible  23  106  175  175  63  60  65  Melnica  22  109  173  173  63  60  71  Toplek 2  21  174  175  175  88  73  64 considering the thermally altered parts found in the nearby dolomite marble formation north from Melnica ( Fig. 2; Temovski 2016). At Karši Podot the Mg/Ca ratio is similar, although the water is discharging from dolomite marble. This could be due to incongruent dissolution within the dolomite formation, as such a process has been identified in Karši Podot Cave. Within the dolomite marble, the porosity was formed by ghost-rock weathering (e.g. Dubois et al. 2014), i.e. preferential dissolution of calcite due to cooling of thermal waters, leaving in-situ dolomitic sand residue (Temovski 2016(Temovski , 2017. As the slowly moving thermal waters were lacking in sufficient energy to carry away the remaining dolomite sand, the penetrable-size cave passages were formed only when the high-energy back-flooding water of Crna Reka removed the dolomite sand residue. The higher Mg/Ca ratio at Gugjakovo indicates that the source of the fresh water component is the dolomite marble formation to the north, instead of the previously considered Upper Cretaceous limestone formation (Temovski 2016).

The groundwater flow systems
Compared to what is commonly reported for thermal karst groundwater at the discharge areas (Goldscheider et al. 2010), the studied springs generally have lower total dissolved content (<1,000 mg/L) and low values for dissolved SO 4 2− and Cl − . Lukewarm spring waters (<30°C) of lower salinity discharging at the Buda Thermal Karst (Mádl-Szőnyi and Tóth 2015) have been interpreted as part of an intermediate groundwater flow in unconfined settings, that has somewhat increased Cl − concentration due to mixing along the discharge zone with basinal fluids of a confined, deeper, regional groundwater flow system. Groundwater of lower salinity related to hypogene karst systems is also found in regions where young volcanism is present, generally associated with hydrothermal activity, with higher concentrations of CO 2 and/or H 2 S (Klimchouk 2017). For the Mariovo hypogene karst systems, based on the geochemical composition of the studied springs, mixing of waters from two groundwater flow systems is identified: a young cold freshwater and an older, chemically more evolved thermal water. The mean residence time of the young shallow groundwater is constrained by the 3 H-3 He apparent ages, and ranges from recent (<3 years) at Toplek 4,~10 years at Manastir, 30 years at Podot locality and~50 years at Toplek 2. Manastir and Toplek 4 have the youngest groundwater, which reflects the relatively fast groundwater flow of the local groundwater systems developed in highly permeable travertine and dolomite of higher permeability, respectively. The young groundwater at Podot and Toplek is somewhat older, and reflects either slower groundwater flow of local groundwater systems in less permeable marble and dolomite rocks, or possibly (e.g. in the north stretching marble stipe at Podot) an intermediate groundwater flow system.
The 14 C age of the old groundwater component at Podot, if based on the calculated carbon isotope composition (~30 ka), is clearly overestimated. This is due to addition of 14 C-free carbon, partly due to dissolution of carbonate rock, but more importantly due to the large contribution (up to 54% of total DIC) of endogenic CO 2 (Table 6). Nevertheless, the deadcarbon dilution-corrected 14 C age at Podot (~15 ka) still reflects long mean residence time of the old groundwater, with the estimated mean residence time at Melnica between~6 and 15 ka. The circulation depth of 0.9-1.2 km at Melnica and Podot indicates deeper circulation of the old groundwater, which can be achieved if the flow is along the dip of the carbonate formation. Considering the topographic position of the discharge areas, the local geological constrains, the estimated circulation depth, the inferred water-rock interactions and the higher recharge elevations, the old groundwater component at Podot represents the discharge of a regional groundwater flow system. The groundwater at Melnica has somewhat lower residence time and represents either a shorter flow path (a subsystem) of the same deeper regional groundwater flow system after reaching a deep fault, or a shallower intermediate groundwater flow system. However, the latter is probably less likely, considering the very similar geochemical properties with the old groundwater at Podot. The very high uncertainty in the 14 C age of the old groundwater at Toplek does not allow for any meaningful interpretation.

The hypogene speleogenetic settings
Most of the carbonate rock dissolution likely occurs along the rising limb of the flow in the discharge zone. While increased dissolution here due to mixing with shallower groundwater (i.e. mixing corrosion) is possible, the high dissolved CO 2 concentrations are probably of higher importance (Dublyansky 2000). Additionally, cooling due to rising of the thermal groundwater and/or mixing with the cold shallow water, can further increase dissolution due to the retrograde calcite solubility (Andre and Rajaram 2005). These processes are consistent with the findings from Provalata Cave, where the older speleogenetic phase is associated with dissolution of the marble bedrock due to cooling of CO 2 -rich thermal waters (Temovski et al. 2013). Calcite deposited in isotopic equilibrium with the groundwater DIC at Melnica Spring will have a δ 13 C of ca. +4‰, which is comparable with the measured high δ 13 C values (ca. +7‰) in the calcite deposits of Provalata Cave (Temovski et al. 2013). This suggests that this system has operated for a long time, and the abundant travertine deposits in the surroundings (Temovski 2016), deposited in different environments (e.g. lacustrine limestone, paludal limestone, spring tufa), are also likely related to increased alkalinity supplied by these deep groundwater systems.
Mantle helium contribution in the dissolved gases at Podot and Melnica of up to 10% is consistent with the findings on the southern foothill side of Kožuf-Kozjak volcanic system in Greece, where up to 16% mantle helium was reported (Daskalopoulou et al. 2018). H 2 S oxidation is a likely source of the dissolved sulfate at Melnica and Karši Podot. A magmatic origin of H 2 S is also possible based on the sulfur isotopic composition, as well as the presence of mantle helium. Mariovo hypogene karst system shows some geochemical similarities to other hypogene karst systems related to young volcanism such as the presence of mantle derived gases, e.g. Sistema Zacatón in Mexico (Gary and Sharp 2006), Mt. Gambier in Australia (Webb et al. 2010) or Konya Basin in Turkey (Bayari et al. 2009b), except that the CO 2 at Mariovo is of dominantly metamorphic origin. However, it shares none of the morphological expressions of the discharge zone, with large collapse dolines developed at Sistema Zacatón, Mt. Gambier and Konya Basin, and only relatively small caves found in Mariovo. One reason for this can be the structural difference, with Mariovo hypogene karst developed in highly dipping strata, while the others are developed in subhorizontal strata, but also the high difference in fluid flux, with Mariovo having both smaller flow rate and gas concentrations.

Conceptual model of the Mariovo hypogene karst system
The results from this geochemical study allow one to put some constraints on the Mariovo hypogene karst system and develop a conceptual model (Fig. 15). The water is of meteoric origin, likely infiltrated along fault structures in the southern parts of the marble stripe. It reached depth of around 1 km, which could be achieved by circulating along fault structures and/or along the dip of the carbonate formation, and then following northward in strike direction. The output zones are at the intersection of low topography and major fault structures at Podot and Melnica localities. Podot locality, is the furthest (~25 km) and lowest output, with water emerging~15 ka after recharge. Discharge of endogenic gases, mostly metamorphic CO 2 and especially some contribution of mantle helium in the output zone indicates deep setting of these fault structures, related to the extensional tectonics of the area. Melnica is likely a subsystem of the same groundwater system, with a shorter flow path (~18 km) and younger age, formed because part of the groundwater was likely forced upward after reaching a deep fault with a higher gas flux.
Most of the carbonate rock dissolution is likely achieved along the rising limb, due to cooling with acidity acquired from endogenic CO 2 . At depth, closer to the recharge area, the groundwater has contact with rocks from the Kožuf-Kozjak Neogene-Quaternary volcanism. Contact with the underlying metamorphic basement and the granitoid bodies perching it, is also achieved along the northward flow either at depth or along the rising limb in the output zone. In the southern parts, between the recharge zone and Melnica, a schist formation is separating the lower, mostly dolomite marble formation, from the upper, calcite marble formation, partly confining the groundwater system, and forcing the northward deep flow. At the output zones, the deep circulating old water of the regional to intermediate hypogene groundwater system is mixing with the shallow recent water part of the local epigene groundwater systems. This hypogene groundwater system has been active for a long time, as evidenced by the >1.5 Ma age of deposits in Provalata Cave (Temovski et al. 2013). The large travertine deposits found within, or topping, the Pliocene-Quaternary basin sediments near Melnica or Podot are likely related to older stages of the same hypogene karst system.

Conclusion
This study presents an approach whereby various geochemical methods are applied in order to identify different components Fig. 15 Conceptual model of the Mariovo hypogene groundwater system. The NNW-SSE direction of the cross-section represents the strike of the carbonate formation. The thickness of the carbonate formation is the apparent thickness along the dip of the formation and properties of groundwater associated with hypogene karst systems. It demonstrates the use of radiogenic and stable isotopes, noble gases and major ion and trace element composition in a case where sampling was restricted only to the spring sites, and in a limited number of locations. In the absence of clearly defined end-member compositions, by combining multiple geochemical methods the groundwater components are identified, as well as their composition and mean residence time, and their geochemical evolution. Furthermore, the carbon sources contributing to the DIC are identified, based on which the carbon isotope composition of the endogenic CO 2 is modeled. This is then used to correct the radiocarbon composition for the dilution of the soil 14 C signature by addition of 14 C-free carbon from dissolution of carbonate rock and endogenic CO 2 , and calculate a radiocarbon age for the old component.
The obtained results show that the main studied springs represent an output part of a regional hypogene karst groundwater system with a deep-circulating (~1 km), old (~15 ka), thermal (≥60°C) water, where the old groundwater mixes with young (<50 years), cold (<14°C) and shallow epigene karst groundwater. These output zones of the hypogene karst system are located at the interception of low topography and deep fault structures along which there is circulation of deepseated gases, dominantly CO 2 of metamorphic origin (+4.5‰ δ 13 C), with a contribution of mantle helium. The chemical composition of the water indicates interaction at depth with volcanic rocks from the N-Q Kožuf-Kozjak volcanism, as well as with metamorphic rocks and granitoids of the Pelagonian basement. The groundwater geochemistry results of this study complement previous findings based on geochemistry of fossil cave deposits, and confirm that dissolution of carbonate rocks due to cooling of CO 2 -rich thermal waters is the main hypogene geochemical process in the area.