Identification of Natural and Anthropogenic Geochemical Processes Determining the Groundwater Quality in Port del Comte High Mountain Karst Aquifer (SE, Pyrenees)

The Port del Comte Massif (SE, Pyrenees) contains one of the most important vulnerable and strategic karst aquifers for supplying freshwater to the city of Barcelona (Spain). It is a fragile system, whose possible environmental impact is highly conditioned by land use. To improve the hydrogeological knowledge of the system, between September 2013 and October 2015, a detailed fieldwork was carried out for the revision of the geological model, the inventory of water points, and the in situ physico-chemical characterization on major elements and isotopes of up to a total of 43 springs, as well as precipitation water. This paper focuses on the characterization of the geochemical processes that allow explanation of the observed chemical variability of groundwater drained by the pristine aquifer system to determine the origin of salinity. The results show that the main process is the dissolution of calcite and dolomite, followed by gypsum and halite, and a minor cation exchange-like process. Sulfur and oxygen isotopes from dissolved sulfate in the studied springs point out a geogenic origin related to the dissolution of gypsum from Triassic and Tertiary materials, and that the contribution from anthropogenic sources, like fertilizers, is lower. Nitrate in groundwater is not an important issue, with a few localized cases related with agricultural activities. The multidisciplinary approach has allowed the development of a consistent hydrogeological conceptual model of the functioning of the aquifer system, which can be replicated in other places to understand the geogenic character of the hydrogeochemistry.


Introduction
High mountain karst aquifers are strategic freshwater reservoirs to maintain dependent ecosystems downstream, many of which are in semiarid zones. Globally, 68.9% of all the surface exposures of carbonate rocks occur in hills and mountain areas. Roughly, 20-25% of the world's population depends directly or indirectly on water supplies from karst aquifers [1]. Given the relevance, it is essential to characterize such mountain karst aquifers and protect them to avoid undesirable quality issues in the stored water resources, because karst aquifers have been shown to be very vulnerable, especially to pollution, given and limited self-depuration capacity of such aquifers make them especially vulnerable to pollution [2].
The aim of this work is to point out the geogenic origin of solutes as well as possible anthropogenic contributions driving the hydrogeochemical composition of GW in a pristine rural karst system to form a comprehensive picture of the hydrogeochemical behavior of this vulnerable type of aquifer in order to prevent the effect of changes due to future scenarios for global change.

Geographical and Climatological Settings
The PCM is approximately 100 km north of Barcelona (Spain), in the south of the eastern Pyrenees ( Figure 1). The elevation of the mountainous massif ranges from 900 m a.s.l. to 2383 m a.s.l. at the 'Padró dels Quatre Batlles' peak. The study area covers an extension of about 110 km 2 . The NW and SW part of the massif drain to the Segre River basin, whereas the eastern and southern parts drain to the Cardener river. The main source of this last river is the Fonts del Cardener spring (M-22, Figure 1). The area is covered by forest (63.9%) and mountain meadows (21.1%), which include one alpine ski resort and one cross-country (Nordic) ski resort, areas without soil and vegetation (11.5%), agricultural cultivation areas (1.5%), and the rest corresponding to residential zones (0.5%) [29]. From a climatological point of view, and according to the Köppen-Geiger classification system [40], the study area is characterized by a cold climate without a dry season and with a temperate summer. Within the study area, there are two meteorological stations located at 2315 m a.s.l. (SMC) and 1800 m a.s.l. (AEMET). For the period 2005 to 2019, the average annual precipitation, temperature, and potential evapotranspiration (Hargreaves) measured at the SMC station are 1055 mm, 3.24 • C, and 525 mm, respectively. During 3 to 4 months in the winter, depending on the year, the snow cover is above 1800 m a.s.l. The infiltration of the snowmelt during the spring season generates the main contribution to aquifer recharge in the PCM [29]. The massif presents a characteristic karst pattern in the upper part, which appears in the form of sinkholes and karren fields that allows a fast infiltration of the meteoric waters and prevents large surface runoff events.

Geology and Hydrogeology Setting
From the geological point of view, the PCM constitutes an independent thrust sheet limited by other thrusts with complex structural relationships ( Figure 1). Internally, the massif presents a set of folds (anticlines and synclines) and some faults. The folds have a characteristic NE-SW direction parallel to the NW boundary of the thrust sheet [31]. Stratigraphically, the massif contains materials from the Triassic (Muschelkalk limestones and Keuper evaporites), the Cretaceous (limestones, calcarenites, shales), and the Paleogene from the Eocene to the Oligocene (karstified limestones with dolostones, sandstones, and marls). The Jurassic carbonates only outcrops in the NW part of the study area. The massif has a total thickness exceeding 1300 m. From the geomorphological perspective, the PCM is characterized by presenting a rounded or flat relief in the highest part, where no vegetation cover is present and almost without any type of soil development exists. The rest of the massif is covered by mountain meadows and forests, with little thin soil. The area most affected by karstification becomes visible progressively, from 1950 m a.s.l. to the top. shales, marls and limestone (Kgp), multicoloured 'redbed' facies clay deposits; (6) Lower Eocene: fissured/karstified alveoline limestones and dolostones (PPEc), and includes colluvial Quaternary formations that partially overlap (Qpe, Qvl); (7) Lower Eocene: marls, sandstones, and limestones (PEci); (8) Lower Eocene: fissured/karstified micritic and bioclastic limestones (PEcp1, PEcp2); (9) Middle Eocene: sandstones, marls, conglomerates, limestones, and evaporites (PEalb, PEm1, PEmb and PExb), including colluvial Quaternary deposits (Qcoo) and alluvium (Qoo) that partially overlap; (10) Upper Eocene: alluvial systems: conglomerates and sandstones; and (11) Oligocene alluvial system: conglomerates and breccias deposits and sandstones (POcgs, POmlg). The breccia deposits covering the Lower Eocene in the upper part of the massif are very thin. The epigraphs in parentheses correspond to the geological units [41] where the springs included in this work are located (modified from Herms et al., 2021).
Hydrogeologically, the PCM can be considered an independent unit and a multiaquifer system. The Lower Eocene fissured and karstified limestones and dolostones form the main aquifer. It constitutes one of the most important high mountain karst aquifers (HMKA) of the Catalan Pyrenees. The other aquifers and aquitards in the study zone are related to the Cretaceous limestone formations, Triassic limestone and evaporites, other Paleogene conglomerates and sandstones, and small local quaternary aquifers (draining small areas) located at low or medium elevations. The Garumnian materials (Upper Cretaceous-Lower Paleogene), composed of shales, marls, and multicolored clay deposits (geological unit 5 in Figure 1), constitute a low permeability layer that acts as a lower impervious limit for both the Tertiary materials and the main karst aquifer. In addition, the geological structure of the system strongly influences the location and hydrodynamical behavior of the springs and their geochemical fingerprint.
Although many high mountain karst systems play a strategic role in terms of availability of underground water resources, such as the case of PCM [29], frequently they are not sufficiently known because the conventional hydrogeological investigation techniques are often difficult to be applied given the complicated access and harsh working conditions typically existing in high mountain zones [42][43][44][45][46].
From a hydrological point of view, there are relevant processes, such as the aquifer recharge or the mean transit time, that need to be characterized to correctly manage the water resources generated in such alpine groundwater systems. Spring hydrograph analysis and the use of environmental tracer methods allow the characterization of aquifer recharge and discharge processes, assessing spring vulnerability, as well as estimating the available water resources [47][48][49][50]. Rainfall-runoff hydrological models are useful to simulate the behavior of such complex mountain karst systems. Such models can be broadly categorized into lumped, semi-distributed, and fully-distributed models [51]. Semidistributed and lumped parameter models (LPMs) are often used to simulate the behavior of high mountain aquifer systems, because they do not require a detailed hydrological knowledge of the physical system, and therefore, they are specially well suited when the hydrogeological systems are poorly characterized. Additionally, the stable isotopes of water (δ 18 O and δ 2 H) in precipitation have proved to be good environmental tracers for investigating the dynamics of such hydrological systems karst systems [52]. These tracers enter the system as recharge, migrate downslope exploring the entire hydrological system, and leave the karst aquifer through the springs discharge. In the case of PCM, [29] presents a hydrogeological characterization of the aquifer system response that includes: (1) The hydrodynamic behavior of the system, simulating the system response with a set of semi-distributed rainfall-runoff HBV models [53,54], while taking into account the elevation dependences of both the hydrometeorological variables (e.g., precipitation and temperature) and the related processes (e.g., snow accumulation and ablation). The estimated groundwater storage capacity of the system is 35.2 hm 3 , and the mean annual groundwater discharge is 15.4 hm 3 ; and (2) The estimation of the mean transit times corresponding to the main springs draining the aquifer system. This is done by using a set of LPMs models [55] to simulate the environmental tracers' content evolution in groundwater. The LPMs were implemented for the most important karst springs of the PCM systems, i.e., the four 'regional springs' named as M-22, M-25, M-31, and M-43 in Figure 1. The results indicate that the PCM karst system presents a relatively short mean transit time (~2.25 year). This result is relevant if the hydrological high conductive features existing in the karst system are taken into account, which may favor a fast contaminant migration from the recharge to the discharge areas in the case of eventual surface spills of contaminants.
In the framework of the 'GeoERA Resources of groundwater harmonized at cross-border and pan-European scale (RESOURCE) Project' (2018-2021) funded from the European Union's Horizon 2020 research and innovation programme, the PCM karst aquifer was incorporated as a pilot case study among others across Europe [56]. In this framework, a well-known hydrodynamic typology classification method proposed by Mangin (1975) is applied [57], which is based on a recession curve analysis of the spring hydrograph. The method describes the hydrodynamic behavior of a karst system as a function of two indices, k and i, defining the extent of the karst phreatic zone and characterizing the infiltration conditions, respectively. In the case of spring M-22, and for a diagram k vs. i [57], the indices lie in the k < 0.5 and i > 0.5 domain, thus indicating a complex, large karst system, made up of several sub-systems. This result is consistent with those obtained for other important karst springs in the Pyrenees, such as Fonts del Llobregat springs [58], which are one of the main groundwater resource for the Barcelona metropolitan area. This manuscript presents the work done to characterize the PCM system from the hydrogeochemical point of view. In this sense, several geochemical fieldworks were carried out on 43 springs during the period September 2013-October 2015. The 43 springs were grouped from a geochemical point of view [59] applying an approach based on a the Gaussian mixture model (GMM) [60], obtaining four groups ( Figure 2). GMM is a 'soft' model-based clustering method that avoids the degree of subjectivity assumed by the classical 'hard'-or heuristic-based-algorithms (e.g., k-means, hierarchical clustering) [61] to determine the relevant number of clusters. Following [59], the four spring clusters found in PCM can be summarized as:

•
Cluster A: 27 springs characterized by low mineralization and dominated by slightly alkaline Ca-HCO 3 water type, which is associated with the Eocene carbonate materials conforming the main aquifer of PCM. • Cluster B: 10 springs that include different types of water from Ca-HCO 3 to Ca-HCO 3 -SO 4 , Ca-SO 4 -HCO 3 , and Ca-SO 4 , which are characterized by moderate mineralization. These springs are located both inside and outside the structural limits of the PCM trust sheet. The springs located inside the limits are mainly found in materials from the Cretaceous and Triassic (Keuper) that outcrop in the area. These materials underlie the main aquifer of the massif (the Eocene karst carbonate system). In the southeastern part of the study zone, there are five springs related to sediments with a high content of tertiary gypsum from the Eocene-Oligocene Beuda's gypsum Formation, pinched out within the South Pyrenees thrust fault in the front SE of the PCM. • Cluster C: 4 springs with water types of Ca-HCO 3 and Ca-HCO 3 -Cl. Three of these springs are located at the boundaries of the PCM sheet. • Cluster D: corresponds to two salty springs with Na-Cl facies that are in the eastern and western limits of the PCM thrust sheet, respectively. They are characterized by remarkably high mineralization and are saturated relative to gypsum. According to Herms et al. (2019) [29], the recharge of the main aquifer of the PCM is produced by diffuse infiltration of precipitation-rainfall and snowmelt, but it is also concentrated through the well-developed karstic elements, mostly situated at the top of the massif. The infiltrated water flows vertically through the unsaturated zone, which may be thicker than 1000 m, towards the saturated zone of the aquifer. GW discharges through a large network of springs. In the framework of this work, more than 100 springs have been identified in the study zone. Most of them discharge local sub-surface waters with low flow rates ranging between <0.1 L/s and 10 L/s. There are four 'regional' springs (M-22, M-25, M-31, and M-43; Figure 1) with mean flow rates from 1 L/s to 900 L/s, whose recharge areas are at medium to high elevation areas. Springs M-22 and M-25 discharge through Quaternary deposits overlying the Lower Eocene limestones, M-31 discharge directly through the limestone outcrops, and M-43 through well-developed karstic conduits affecting the conglomeratic materials of the Ebro Basin (the southern foreland basin of the Pyrenees), which are located just at the southern border of the PCM mantle. There are also some GW diffuse discharge zones, especially in the northern sector of the PCM. The regional water table main aquifer is located at elevations between 1000 and 1100 m a.s.l. [29].

Field Measurements, Sampling, and Laboratory Analysis
A total of 43 springs (Figure 2A) were monitored in this research for the period from September 2013-October 2015, in which the springs were sampled twice per year (i.e., before snowfall in October and after snowmelt in April). Additionally, springs M-04, M-20, M-22, M-25, M-31, and M-43 were regularly sampled more frequently, every three to four weeks, along the same sampling period. In every case, the "in situ" physicochemical parameters (pH, EC, T, redox, alkalinity, TDS) were measured. A total of 288 GW samples were obtained. Table A1 in Appendix A summarizes the details of the GW sampling campaigns conducted in this work. Tables A1-A3 in Appendix A provide the chemical characteristics of major constituents for the 43 water springs (median values for the whole campaigns carried out from September 2013-October 2015) and the rest of samples. Hereinafter, for simplicity, the different ions are indicated without the charge if this does not create confusion.
The hydrochemical composition of precipitation is obtained from an open (bulk) precipitation gauge [62] at the neighboring meteorological station of la Molina (42 • 20 30" N, 1 • 57 14" E, altitude 1704 m a.s.l.), which is located 30 km to the NW of PCM. The hydrochemical composition of recharge ( Figure 2B) is estimated from that of precipitation but applying an evapo-concentration process to simulate the effect of actual evapotranspiration in the PCM as estimated by [29].
To estimate the impact of evapo-concentration in the meteoric water percolating through the soil, the modelling results obtained by Herms et al. (2019) [29] with the HBV model [63] for the springs M-25, M-43, M-31, M-22, and M-20 are used (springs in Figure 2). A regression line between the evaporation factor (Ef) and the recharge zone elevation of these 5 springs is obtained, thus allowing estimation of the associated Ef for the other 43 springs of the study zone. The average calculated value for the recharge chloride content is 2.2 mg/L, which is consistent with that obtained with the chloride mass balance method by [64]. Figure 2B and Table A5 in Appendix A show the hydrochemical composition of precipitation [62] and the estimated average recharge evapo-concentrated water chemistry in the PMC applying a reduced concentration factor as estimated by Herms et al. (2019) [29].
To characterize the temporal variation of the isotopic content of precipitation at the PCM, precipitation was sampled seasonally in 8 cumulative rain gauges installed at elevations between 896 and 1935 m a.s.l. A total of 71 precipitation water samples were taken, besides 10 snow samples (3 of them corresponding to artificial snow samples from a ski resort situated within the catchment area), and 2 surface water samples, from water ponds used for manufacturing artificial snow. In addition, rock samples containing gypsum were taken to characterize isotopically local sulfate. Figure A1 and Table A7 in Appendix A show the location of the sampling points and the results obtained.
All water samples were filtered in the field using a 0.45 µm membrane filter and stored in new 200-500 mL polyethylene bottles washed with diluted nitric acid and rinsed with the water to be sampled prior to sampling. Samples for cation analysis were acidified to pH<2 with ultrapure HNO 3 to prevent precipitation, while samples for anion analysis were not acidified. Samples were preserved at 4 • C until laboratory measurement. The physico-chemical parameters of GW (T, CE, pH, Eh, and TDS) were measured in situ by a portable Hanna meter (Multiparameter Water Quality Meter HI9829). Total alkalinity was determined in situ by using the titration method in the first four campaigns as well as with an Alkalinity Test Checker de (HI755) of Hanna Instruments.
The major ions, cations (Ca , Mg, Na, K, NH 4 ), and anions (Cl, NO 3 , HCO 3 , CO 3 , SO 4 were measured in the Laboratori Ambiental d'Aigües de Terrassa: the cations were analyzed by inductively coupled plasma atomic emission spectrometry-ICP-OES (Agilent 5100 DV), except ammonium that was determined by ultraviolet-visible (UV_VIS) spectrophotometer, and the anions by ion chromatography (Dionex, DX-120). Table A2 Appendix A summarizes the median concentrations for the 8 major ions of the total 43 monitored springs.
The isotopic composition (δ 2 H and δ 18 O) of the low-salinity water samples was determined in the Center of Hydrogeology of the University of Málaga (CEHIUMA), with a Picarro ® "L2130-I" cavity ring down-spectroscopy analyzer, which is based on cavityenhanced, near-infrared laser absorption spectroscopy procedures, tuned on a narrow spectral region. The analytical uncertainties for δ 2 H and δ 18 O are ±0.2 ‰ and ±1.0 ‰, respectively. According to Coplen (2011) [65], several international and laboratory standards have been interspersed for normalization of analyses. The standards used (WICO-13, WICO-14, WICO-15) were calibrated in an interlaboratory comparison exercise [66]. In the case of high-salinity water samples, analysis of δ 2 H and δ 18 [69]. Solid gypsum samples from Triassic and Tertiary evaporites were previously dissolved in water. The δ 34 S was analyzed in a Carlo Erba Elemental Analyser (EA) coupled in continuous flow to a Finnigan Delta XP IRMS. The δ 18 O was analyzed in duplicate with a Thermo Quest high-temperature conversion analyser (TC/EA) unit in continuous flow to a Finnigan Matt Delta XP IRMS.
To check the accuracy of the analytical results, ionic balance errors were calculated using the USGS software PHREEQC ® [70] using the phreeqc.dat database for all water springs, except for the brines from M-30 and M-41. Most of the samples have ionic balance errors below the recommended standard of ±5% [71]. Isotope ratios were calculated using both international and internal laboratory standards. Notation was expressed in terms of delta (δ)‰ relative to the international standards: V-SMOW for δ 18 O and δ 2 H, atmospheric N 2 for δ 15 N, and V-CDT for δ 34 S. The precision of the analyses calculated from the reproducibility of standards interspersed in the analytical batches was ±0.3‰ for δ 15 N, ±0.2‰ for δ 34 S, and ±0.5‰ for δ 18 O of SO 4 and NO 3 .

Application of the Dual-Isotope Approach for δ 34 S and δ 15 N
The existence of NO 3 and SO 4 in GW may pose an environmental risk in many mountains and rural areas with pristine waters. The use of stable isotopes of N and O of dissolved NO 3 and S and O of dissolved SO 4 , together with the geochemical data, has proven to be a useful tool to evaluate the origin of solutes [7,20,72,73]. These tools help to improve the knowledge of the hydrogeological system, but also to understand both the natural and anthropogenic geochemical processes driving the GW quality in the aquifer system. In this research, stable isotopes of δ 15 N NO3 and δ 18 O NO3 as well as δ 34 S SO4 and δ 18 O SO4 were used to identify NO 3 and SO 4 sources in groundwater using the well-known dual isotope approach [74]. Ranges of NO 3 and SO 4 isotope compositions of the main potential sources were obtained from [73].

Determination of Proportional Contributions of NO 3 and SO 4 Sources
To estimate the relative contribution of different sources, stable isotope mixing models have become a common tool in environmental studies. Beyond the well-known limitations of the classical mass-balance mixing models-related to the restriction of taking at maximum n + 1 sources for n isotopes [75]-nowadays, the Bayesian isotope mixing models (BMMs) [76,77] are the focus of attention to determine the probable source apportionment. BMMs are traditionally used to identify biogeochemical sources. They allow estimation of the probability distribution for the relative contribution of each source considering the uncertainty associated within the sources themselves and their isotopic compositions. Since a few years ago, these techniques have been used to assess the contribution of NO 3 pollution in general studies [78][79][80][81][82][83][84]. They have also been used specifically in karst studies [85,86], and in very few cases for SO 4 studies [87,88]. Different BMM codes for R and MATLAB have been reported, among them: SIAR, MixSIR, SIMMR, and MixSIAR. In this research, the Stable Isotope Mixing Model (SIMMR) package for R, an updated version of the Bayesian isotope mixing model named as SIAR [76,77], was used to determine the proportional contributions of natural and anthropogenic NO 3 and SO 4 sources into groundwater. The SIAR model (Stable Isotope Analysis in R) is an open-source software package for R. It uses a Markov chain Monte Carlo method to simulate plausible source proportions. The SIAR model is formulated according to Equations (1)-(4): where X ij is the isotope value j of the mixture i, in which i = 1, 2, 3, . . . , N and j = 1, 2, 3, . . . , J; S jk is the source value k on the isotope j (k = 1, 2, 3, . . . , K) and is normally distributed with mean µ jk and standard deviation ω 2 kj ; p k is the proportion of source k, which is estimated by the SIAR model; c jk is the fractionation factor for isotope j on source k and is normally distributed with mean λ jk and standard deviation τ 2 kj ; and finally, ε ij is the residual error representing the additional unquantified variation between individual mixtures and is assumed to be normally distributed with mean 0 and standard deviation σ 2 j .

Delineation of the Main Recharge-Discharge Pathways
The stable isotopes of water δ 2 H H20 and δ 18 O H20 values were used by [29] to study the response of the hydrologic system to the seasonal variation of the isotope content in the recharge waters, estimating the Local Meteoric Water Line (LWML) and the local Isotopic Altitudinal Lines (IALs) for δ 2 H H20 and δ 18 O H20 . In this work, the IAL is used to support the building of the hydrogeochemical conceptual model. To this end, the main Recharge-Discharge Pathways (RDP) of the system are delineated based on (1) the estimated centroid of the recharge zone elevation associated with every spring, and (2) the geological structure of the PCM.

Inverse Hydrogeochemical Modeling for the Quantification of Chemical Processes
Using the PHREEQC code [70], mass-balance inverse geochemical models are applied to analyze the chemical changes that occur from the recharging areas to the discharge points, and to validate the conceptual hydrogeochemical model of the PCM. Inverse modelling is based on a geochemical mole-balance model, which calculates the transfer of moles of minerals and gases that must enter or leave a solution, while accounting for the differences in an initial and a final water composition along a hypothetic GW flow line. In this work, the study is focused on 14 springs that are considered representative of the main Recharge-Discharge Pathways (RDP) of the system, which were previously defined based on the recharge elevations inferred by means of the IALs. The differences in hydrogeochemical composition between the springs are assumed to be due exclusively to reactions between GW and the minerals within the PCM. The selection of the solid phase reactants is based on the geological knowledge of the main lithologies, which comprise calcite, dolomite, gypsum, and halite. In addition, soil gas CO 2 is also considered. The existence of marl and clay materials in the limits of the PCM, as in the Cretaceous and Triassic materials, suggests that ionic exchange-like processes between cations Ca, Mg, Na, and K and an exchanger X might occur. As all water samples are undersaturated according to the calculated saturation index, the inverse modeling was constrained so that dolomite, gypsum, and halite only dissolve whereas calcite is allowed to both dissolve and precipitate. In summary, the expected reactions responsible for the groundwater composition can be defined as:

•
Dissolution of carbonate minerals, such as calcite and dolomite, and precipitation of calcite according to Equations (5) and (6) [89]. Equation (7) is obtained as the sum of Equations (5) and (6), and it shows that the molar ratio Ca 2+ /Mg 2+ is 3:1: • Dissolution of evaporite minerals, such as gypsum and halite, according to Equations (8) and (9): • Dedolomitization processes according to Equation (10) [71], which causes an increment of Ca 2+ due to gypsum dissolution (as indicated in Equation (8)) and precipitation of calcite: • Ion exchange reactions due to weathering reactions in marls, shales, and clays associated with Triassic and Cretaceous layers, according to the following Equations (11)- (14): Silicate minerals in the aquifer were generated in the marine environment in which the carbonates were deposited, and this affects the sorbed cation composition. Dispersed silicates in the carbonate rock matrix may progressively contribute Na + as carbonates are dissolved, as well as a fraction of the K + , at the time that HCO 3 increases. This is a minor process. Marl and clay layers may contain evaporitic minerals that slowly diffuse out, contributing Cl − , SO 4 2− , and cations that correspond approximately to halite and gypsum dissolution, although with a modified cation composition. All these processes can be assumed steady in a hydrogeological system under natural conditions, so cation exchange is minor and not significant. However, the silicate weathering and porewater diffusion in clays produce similar results as cation exchange, and therefore they are considered as such for general treatment and called cation exchange-like process.
There are no data on soil CO 2 partial pressure and therefore this value and the carbon isotopic composition must be assumed. Then, the CO 2 dissolved in the meteoric water recharging the aquifer is assumed to be (1) in equilibrium with atmospheric CO 2 partial pressure for elevations above 1900 m a.s.l. and equal to log(PCO 2 ) = −3.2 (no vegetation), and (2) equilibrated with the soil CO 2 for elevations below 1900 m a.s.l., with the CO 2 content estimated by the following equation [71] to consider the existence of decaying vegetation and root respiration: where PCO 2 corresponds to the mean growing-season soil CO 2 partial pressure and ETP is the mean potential evapotranspiration. The mean ETP in the PCM is 525 mm/year [90] and the obtained log(PCO 2 ) is −2.23.
Considering the geological structure of the PCM along with the estimated average recharge altitude associated to every sampled spring, a total of 14 representative Recharge-Discharge Pathways, named as RDPs 1 to 14 ( Figure 3), were considered for inverse hydrochemical modeling with PHREEQC. Every RDP tries to recreate the GW flow line integrating all the processes driving the hydrochemical composition of the corresponding spring discharge, from which the RDP adopts the cluster type. Therefore, the RDPs are classified in clusters as:

Results and Discussion
The hydrochemical and isotope data of GW corresponding to the different sampling campaigns conducted in this work are reported in Tables A1-A8, Appendix A.

Saturation Indexes
The saturation indexes (SIs) were calculated using PHREEQC ® [70]. The SI relative to calcite for all the dataset ranges from 0 to 0.82, indicating calcite saturation to oversaturation throughout the PCM ( Figure 4A). The SI relative to dolomite ranges from −1.32 to 0.61, except for a 2.  Figure 4C). Figure A4 shows the relationship between SI relative to calcite, dolomite, gypsum, and halite with respect to TDS [mg/L], besides SI of gypsum with respect to SI of halite, and SI of calcite with respect to SI of gypsum. In all cases, it is possible to observe a clear separation between all the clusters A, B, C, and D. In addition, a trend of increasing TDS can be observed from clusters A, C towards cluster B, and later to the extreme D. In general, the springs of cluster A have lower TDS than those of cluster B. Figure 5A shows that despite all the samples being in equilibrium or slightly supersaturated with respect to calcite, cluster A shows much less TDS content than cluster B. Additionally, Figure 5B indicates that most of the samples are subsaturated relative to dolomite, except in four samples from cluster A (springs M-16, M-17, M-18, and M-19 located in the northern part of the PCM that discharge the PPEc unit in contact with the PEci unit), four samples from cluster B (M-28, M-36, M-09, M-13), and one sample from cluster C (M-20) that is close to the equilibrium relative to dolomite. Figure 5C shows how all the springs but one (M-30 from cluster D) are under-saturated with respect to gypsum. Moreover, here the separation between clusters A and C with respect to B and D can be clearly observed, indicating that the samples of cluster B are influenced by Triassic (Keuper) and Tertiary (Beuda Formation) formations that contain gypsum. Furthermore, the clusters' separation is even clearer when looking at the relationship between SI relative to halite and TDS ( Figure 5D), SI of gypsum, and SI of halite ( Figure 5E), and less evident for SI of calcite and SI gypsum ( Figure 5F).

Identification of Hydrogeochemical Processes Explaining the Spring Clusters
To determine the main rock-water interactions within the PCM system driving the hydrogeochemical composition of GW, it is necessary to focus on the relationships between major cations and anions. This will help to decipher the main hydrogeochemical processes conditioning the cluster definition presented by [59]. The data presented in this analysis for every spring corresponds to average content values for the monitored period September 2013-October 2015 (Table A2 Appendix A). The relationship between Ca and SO 4 is depicted in Figure 5A. The samples of cluster B follow the line of slope 1:1, thus indicating most probably the dissolution of gypsum as the origin of sulfate in cluster B, while cluster A and C are not clearly related with this process. This suggests that calcite and dolomite are not the primary sources of Ca for cluster B, but they might be for cluster A and C.
The relative contribution of calcite and dolomite in the carbonate weathering processes can be approached by looking at the molar ratio rCa/rMg (r = concentration in meq/L) ( Figure 5B). A molar ratio of 1 (1:1 slope line) indicates pure dolomite contribution (Equation (6), ratio values between 1 and 3 indicate a dominance of dolomite dissolution with some calcite contribution, rCa/rMg = 3 indicates dissolution of both calcite and dolomite according to Equation (7)). Larger ratio values mean a predominance of calcite dissolution plus certain dolomite contribution, and finally, rCa/rMg = 10 represents a total contribution of calcite, beyond the evident contribution of sulphate dissolution in cluster B, as indicated in Figure 5A. As can be shown, all the GW samples have rCa/rMg > 1. Almost 50% of the molar ratio values range between 3 and 10, and 16% between 1 and 3. Cluster A shows 12 springs with rCa/rMg >10, 11 springs with 3 < rCa/rMg < 10, and four springs with rCa/rMg < 3. The four springs of this latter case (M-16, M-17, M-18, and M-19) also show saturation indices (SI) relative to dolomite > 1. Dolomite is present in the recharge area within the PPEc unit. Cluster B contains two springs (M-10 and M-21) with rCa/rMg > 10, seven springs with 3 < rCa/rMg < 10, and one spring (M-36) with rCa/rMg < 3. Cluster C shows one spring (M-27) with rCa/rMg >10, two springs (M-23 and M-42) with 3 < rCa/rMg < 10, and one spring with rCa/rMg < 3. The two springs of cluster D show rCa/rMg values equal to 0.18 (M-30) and 4.32 (M-41). These results suggest that, except for the four cited samples from cluster A that drain the upper part of the PCM, where the materials associated with the PPEc unit are richer in dolomite, the contribution of dolomite increases as the discharge altitude decreases.
The influence of calcite dissolution on cluster A can be observed looking at the ratio rCa/rHCO 3 in Figure 5C. Here, the samples of the cluster A plot align along the line of slope 1:1, representing the stoichiometric ratio of calcite dissolution. To evaluate the influences of the combined dissolution of calcite and dolomite on karst groundwater chemistry, Figure 5D shows the ratio between rCa + rMg with respect to rHCO 3 . The plot shows that one part of the data for cluster A fits the 3: 1 slope line, suggesting that both dissolution of calcite and dolomite contributes to defining groundwater chemistry, whereas the remainder seems to be concentrated at the base of the 10:1 slope, indicating exclusive contribution from calcite. In the case of cluster C, the four springs plot displaced above while maintaining the slope 1:1, which might be related to ion exchange-like processes adding Ca into dissolution according to Equation (11). The samples of cluster B are clearly scattered relative to the 1:1 slope line, thus confirming that they have other Ca sources than calcite and/or dolomite dissolution.
To confirm that dissolution of carbonates (calcite and dolomite) and evaporites (gypsum and halite) are the dominant processes affecting the hydrochemical features of the different clusters, and that ion exchange is a minor process driving the hydrogeological composition of GW, the relationship between (rCa + rMg) and (rHCO 3 + rSO 4 ) is presented in Figure 5E. Most of the springs match the line 1:1 with just very small shifts for clusters B and C, thus suggesting the existence of a very small ion exchange-like process adding Ca and/or Mg into dissolution (Equations (11) and (12)). Figure 5F shows the scatterplot of rNa content versus rCl. In general, the gravity center of all the clusters is below the line of the 1:1 slope, reflecting a chloride excess that principally comes from the atmospheric chloride deposition, although a minor part, especially in the case of samples from cluster B and C, might be attributed to silicate weathering and ion exchange-like reactions in marls, shales, and clays associated with the Triassic and Cretaceous layers through which the groundwater of these springs interacts.
In this regard, and to get more information about the possible contribution of such an exchange-like process, Schoeller (1965) [91] propose using the chloro-alkaline indices indexes, CAI_1 and CAI_2, of common use in hydrogeochemical studies, which are defined as: Values > 0 of both indexes indicate ion exchange and ion exchange-like processes in which dissolved Na and/or K are retained while Ca and/or Mg are released, or Ca-rich brines are incorporated. Values < 0 point to a reverse ion exchange prevalent process [91] or weathering of alkaline ion-rich silicates. In this line, most of the samples do not show a clear sensitivity of CAI_2 with respect to CAI_1 ( Figure 6A). Nevertheless, the samples of cluster C are clearly located in the first quadrant where both CAI_1 and CAI_2 are >0. Figure 6A shows the plot of (rCa + rMg) − (rSO 4 + rHCO 3 ) versus (rNa + rK − rCl). The dependent variable represents the increment of Ca and Mg that is attributed to processes that exclude weathering by carbonates and evaporites (dissolution of calcite, dolomite, and gypsum) while the independent variable gives information of the increment of Na generated by processes other than halite dissolution. A linear relation between these two variables with a slope equal to −1 indicates the significance of ion exchange-like processes as an important factor controlling the groundwater chemistry and its evolution. In the case of the samples of the PCM, the slope of the regression line between (rCa + rMg) − (rSO 4 + rHCO 3 ) and (rNa + rK − rCl) is −0.90 ( Figure 6B). Such a relationship is basically conditioned by the samples of cluster B and C, which would suggest that there may indeed be ion exchange-like processes between the monovalent ions Na and K, and the bivalent Ca and Mg. Moreover, values of (rCa + rMg) − (rSO 4 + rHCO 3 ) > 0 indicate adsorption of Na and release of Ca [92]. The relationship between both HCO 3 /rNa and rMg/rNa vs. rCa/rNa, in Figure 6C,D, respectively, can be used to check the main geochemical process in the system [15,[93][94][95]. In these figures, it can be shown how the origin of HCO 3 for the samples of cluster A is carbonate dissolution. Cluster C and some samples of Cluster B tend to the silicate weathering zone, which also would be related to sodium-calcium ion exchange-like in the weathering of shales. Besides, the two samples of cluster D, corresponding to the deep flow brines, fit in the evaporite window domain.

Aquifer Recharge Altitude Based on δ 2 H and δ 18 O in Precipitation and GW
The estimation of the recharge elevation associated with the springs sampled during this study is conducted by projecting the mean isotopic content of δ 2 H and δ 18 O associated with every spring discharge to the corresponding IAL obtained by [33], which are shown in Figure 7. Table 1 shows the mean recharge elevation intervals for each cluster. The slope of the isotopic altitudinal line for precipitation (∇ z δ P = ∆δ P /∆z) is −1.9 and −12.1‰/km for δ 18 O and δ 2 H, respectively. The overall isotopic altitudinal line for GW (IAL GW ) shows a slope (∇ z δ GW = ∆δ GW /∆z) of −1.4 and 9.5‰/km for δ 18 O and δ 2 H, respectively. The slope values are larger than those corresponding to the isotopic altitudinal line of precipitation (IAL P ), indicating the existence of aquifer recharge along the mountain slope and mixing at the sampling point, a process also known as slope effect [96].    Figure 7 show the values and graphs of slope of IAL GW for the different clusters. Table A2 in Appendix A summarizes the isotope data for each spring, Figure A3 shows the IAL GW graphs for all water samples, including snow water and water ponds, and Figure A4 shows the location of the sampling points for the 10 snow samples and water ponds. As can be seen in Figure 7, the steepest gradient, and therefore the highest role played by slope recharge, corresponds to cluster A, whose springs show typically a Ca-HCO 3 water composition that is related to the Tertiary karst aquifer, which presents a well-developed epikarst zone, thus favoring the infiltration along the mountain slopes where this Tertiary formation crops out.

Quantification of Hydrogeochemical Processes along the Recharge-Discharge Pathways
The principal results obtained with the application of PHREEQC are depicted in Figure 8, which shows the contribution of the different species to all the 14 RDPs considered in the inverse modelling exercise that was conducted. Table 2 summarizes the complete results for clusters A, B, and C. Additionally, Table A9 Appendix A provides the complete data set of results.    The predominant geochemical process in cluster A is the dissolution of carbonates, mainly calcite (37% of dissolved species), followed by dolomite (9.9%) with a bit part of gypsum (1.1%), halite (2.2%), and a residual part corresponding to ion exchange-like processes (0.8%), which agrees with the observations obtained in the scatterplots of ions. The predominant process in cluster B is dissolution of gypsum (52.3%), followed by dolomite (12.1%), halite (9.9%), and calcite (4%), and a small contribution of ion exchange-like processes (3.2%). In this case, this highlights the great range of gypsum dissolution, with values of 1.64 × 10 −3 mol/L (M-36 in RDP-08 with Ca-SO 4 -HCO 3 water type), 5.23 × 10 −3 mol/L (M-10 in RDP-09 with Ca-SO 4 -HCO 3 water type), 2.11 × 10 −3 mol/L (M-09 in RDP-10 with Ca-SO 4 water type), and 9.94 × 10 −2 mol/L (in the M-28 in RDP-11 with Ca-SO 4 ), which are highly conditioned to the typology and extension of Keuper outcrops with gypsum (Tk) and Tertiary with gypsum (PExb), as well as the transit time of groundwater through these materials. The presence of calcite precipitation with dolomite and gypsum dissolution also stands out, which is an indication of de-dolomitization processes in RDP-08 and RDP-11. This process may be especially important in sample M-28 from RDP-11, which shows SI values of 0.70 with respect to calcite, 0.55 with respect to dolomite, and −0.24 with respect to gypsum [71].
Although it has not been modeled, it is worth noting that, according to Municio (2017) [97], feldspar is encountered in the Cretaceous layers (Kat, KMca units), so a small contribution of feldspar dissolution in the weathering processes might also be possible, to explain the contribution of Na and K [98]. In the case of cluster C, the predominant geochemical process is the dissolution of halite (25.2%), but with a similar contribution of calcite dissolution (24.8%), followed by dolomite (7.6%), with a bit part of gypsum (1.2%) and a reverse ion exchange-like process in the weathering of shales, where the Ca (6.2%) and Mg (3%) ions in the aquifer matrix replace Na in solution that probably comes from halite dissolution.

Identification of SO 4 Sources in GW Based on Stable Isotopes
To infer the origin of SO 4 in GW, the relationship between δ 34 S SO4 and δ 18 O SO4 was considered for the different groundwater samples. In this analysis, the isotopic content obtained for the eight gypsum rock samples collected in the PCM area was included ( Figure 9). The overall mean δ 34 S SO4 and mean δ 18 O SO4 in GW are +6.5‰ and +9.5‰, respectively, with a large variation throughout the study area, between −17.5 and +20.4‰, and between +3.   [105]; and (8) Tertiary evaporites [106]. The long and short dashed red lines define the isotopic fractionation range (ε 34 S/ε 18 O SO4 ) in SO 4 reduction reactions, varying between 2.5 and 4, respectively [107]. Table 3 shows the mean isotopic content by clusters. Tables A6 and A7 Appendix A provide the average values for δ 34   The results suggest that the main source of SO 4 for cluster A might be related to all the factors, from pyrite oxidation to Beuda gypsum dissolution. Due to their intermediate isotopic composition, the role of atmospheric deposition, fertilizers, and sewage cannot be determined. Fertilizers have high sulphate contents [103,108] and cannot be discarded because they are applied in alpine ski resort areas to help the grass on the slopes reach maximum growth, and those applied in agricultural soils, which are planted with potato in some areas of the PCM. The sample with the lowest δ 34 S SO4 content (−17.5‰, in M-38 of Cluster A) is assumed to be, according to Municio (2017) [97], the result of sulfide mineral weathering. The spring M-38 is in the upper limit of the 'sulphide oxidation field' as defined by Van Stempvoort and Krouse (1994) [99]. Here, the sulfide minerals correspond to pyrites from marls with lignite materials, such as those appearing in the The samples of cluster C are mostly related to Triassic evaporites but mixed with other sources. In the case of M-20 (cluster C), the isotopic content is consistent with the existence of a klippe of Triassic with Keuper outcropping, which affects the spring catchment, and with some contributions of fertilizers related to the cross-country (Nordic) ski resort, which is close to the spring. The isotopic composition of springs M-23 and M-42 suggests an origin in the soil. Nevertheless, in the case of M-42, there is an additional contribution of sulfate from atmospheric deposition.
The origin of the isotopic composition of the springs belonging to cluster D is related with the Triassic materials. The deep flow brine in spring M-30 matches perfectly inside the Triassic window, whereas the deep flow brine in spring M-41 falls between the Triassic and the Tertiary window, suggesting that the isotopic composition of this spring is affected by SO 4 contributions from these two origins. This is consistent with the structural geological context in which the spring M-41 is located, where the PCM and the Pedraforca thrust sheets coexist and are related to the ductile materials of the Keuper. As a result, the isotopic fingerprint of the Tertiary gypsum is added to that of the Triassic materials into the GW sampled in this spring.

Proportional Contribution of SO 4 Sources in GW in the PCM
To enhance the existing knowledge regarding the sources of sulfate in groundwater, and to better explain the sources' contribution, a Bayesian isotope mixing model was prepared using the SIMMR package for R (an updated version of the SIAR model [76,77]. Figure 10 presents the corresponding outputs aggregated for the four cluster A, B, C, and D with a horizontal boxplot showing the probabilistic contributions for each SO 4 source. The results of the model indicate that the greatest contribution to springs associated with cluster A, which are mostly recharged in areas with little development of soil cover, comes from fertilizers (proportion~20.4%). They are probably related to their use in (1) the "Port del Comte" alpine ski resort, which is located near the top of the massif (Figure 11), an area that is drained by the four regional springs M-22, M-25, M-31, and M-43, and also other local springs, and (2) in the potato fields that are scattered throughout the massif. The atmospheric deposition also contributes (~18.9%), and finally sulfate from sulphide oxidation (~13.4%). For cluster B, the results confirm the clear effect of geogenic sulfate pollution, normally exceeding the drinking water limits > 250 mg/L of groundwater springs principally located at the lowest parts of the PCM. This is due to dissolution of Tertiary evaporites, mainly in the eastern part, with a mean proportion of~34.4% with respect to the total sulfate contributions, and Triassic evaporites (~29.2%). In cluster C, the model shows a generalized mix of all eight sulfate sources considered, ranging between 10.8% and 13.7%. Regarding cluster D, which is composed of springs M-31 and M-41, the model also confirms that the origin of sulfate is mainly geogenic, related to Triassic and Tertiary evaporites with contributions between 16% and 19.6%. Figure A5, Appendix A, presents the model output for every spring. In relation to cluster A, the results obtained for the spring M-38 suggest that the contribution due to sulfide oxidation (64%) is consistent with the biplot shown in Figure 9B, and it is consistent with the field inspections, where the oxidation of sulphides associated with the Cretaceous limestones located in its drainage basin is observed. Another group of outstanding results from the Bayesian model are those associated with the springs belonging to cluster B, which are located at the eastern end of the PCM, an area where there is Tertiary gypsum outcrop. This is the case of the springs M-21, M-41, and M-40 that present a clear contribution, almost exclusive, from this source, with contributions of 85.2%, 78,1% and 89%, respectively.

Identification of NO 3 Sources in GW and Perspectives on Aquifer Vulnerability in PCM
The most important N cycling reactions in lands are nitrification, mineralization-immobili zation-turnover (MIT), plant uptake, denitrification, and NH 4 volatilization [109][110][111]. MIT refers to the recycling of NO 3 via immobilization as organic N, subsequent mineralization to NH 4 via organic matter degradation, and a turnover back to NO 3 via nitrification [111]. Immobilization, together with plant uptake, are two N assimilation pathways, which involve the production of organic N from inorganic compounds, such as NO 3 , NO 2 , or NH 4 [112].
The GW sampling campaigns revealed the existence of nitrate in some of the springs of the study area ( Figure 11). With a median and a mean value of 3.83 and 7.13 mg/L, respectively, the nitrate content in GW does not seem to be a water quality issue in the PCM. Nevertheless, the spatial distribution of nitrate shows that the sources of nitrate may play a role locally. This is important when recharge is produced in a concentrated way, given that focused recharge facilitates the widespread, rapid incorporation and transport of pollutants to groundwater [3].
From the total 43 sampled springs, only 20 of them showed high enough nitrate concentrations in the GW samples to allow determination of the corresponding isotopic (δ 15 N NO3 and δ 18 O NO3 ) content ( Figure 12A). The isotope content in GW ranges between −0.9‰ and +10.7‰ for δ 15 N NO3 , and between −0.0‰ and +8.2‰ for δ 18 O NO3 , with overall mean isotopic contents of +3.5‰ and +2.6‰ for δ 15 N NO3 and δ 18 O NO3 , respectively. Table 3 shows both the mean nitrate and isotopic content by clusters.  Figure 12A), although some mixing is not discarded. Nevertheless, NO 3 in GW appears to originate from soil organic nitrogen compounds, NH 4 fertilizers, sewage/manure sources, or even from a mixing of them. The highest NO 3 content value is 57.3 mg/L in spring M-32 (cluster A; Figure 11). This spring is located close to a potato crop field, where fertilizers are applied. This is consistent with the results obtained for δ 18 O SO4 and δ 34 S SO4 , which suggest the origin of SO 4 is a mixture of soil and fertilizer sulfate sources contributing to GW (Figures 9 and A5, Appendix A for M-32). There is only one spring, M-28, whose H and O water isotopic composition presents the fingerprint of manure and sewage. In fact, this spring is neighbors a cattle farm where manure stocks are managed.
Concerning the δ 18 O NO3 in the NO 3 from GW, its value depends on the δ 18 O of NO 3 − in water (−11.1‰, SD = 5.08‰) and on the isotopic effect produced during nitrification, which is in turn influenced by the +23.5‰ δ 18 O of dissolved atmospheric O 2 [113] and that of H 2 O. The limited variation of δ 18 O NO3 values along the study seems to indicate a negligible isotopic effect from plant uptake. Additionally, it indicates that denitrification processes were not significant along the studied area.
Most of the samples seem to follow a straight-line relationship between the δ 15 N and δ 18 O, with a factor between 1.3 and 2.1, which is consistent with natural denitrification [114] but with a slight variation. The natural denitrification process would be supported by the negative linear correlation between δ 18 O NO3 and ln (NO 3 /Cl) for the GW samples [115] ( Figure 12B). Nevertheless, in the case of PCM, there is almost no correlation indicating that such a process, if it exists, is not significant.  Figure 11). This spring is located close to a potato crop field, where fertilizers are applied. This is consistent with the results obtained for δ 18 OSO4 and δ 34 SSO4, which suggest the origin of SO4 is a mixture of soil and fertilizer sulfate sources contributing to GW (Figures 9 and A5, Appendix A for M-32). There is only one spring, M-28, whose H and O water isotopic composition presents the fingerprint of manure and sewage. In fact, this spring is neighbors a cattle farm where manure stocks are managed.
Concerning the δ 18 ONO3 in the NO3 from GW, its value depends on the δ 18 O of NO3 − in water (−11.1‰, SD = 5.08‰) and on the isotopic effect produced during nitrification, which is in turn influenced by the +23.5‰ δ 18 O of dissolved atmospheric O2 [113] and that of H2O. The limited variation of δ 18 ONO3 values along the study seems to indicate a negligible isotopic effect from plant uptake. Additionally, it indicates that denitrification processes were not significant along the studied area.
Most of the samples seem to follow a straight-line relationship between the δ 15 N and δ 18 O, with a factor between 1.3 and 2.1, which is consistent with natural denitrification [114] but with a slight variation. The natural denitrification process would be supported by the negative linear correlation between δ 18 ONO3 and ln (NO3/Cl) for the GW samples [115] ( Figure 12B). Nevertheless, in the case of PCM, there is almost no correlation indicating that such a process, if it exists, is not significant. The contamination of groundwater by nitrates in the PCM is mostly related to anthropic activities conducted in aquifer recharge areas. Here, the highest nitrate contents in GW are related to agricultural practices. Other relevant anthropic activities in the study area are restricted to those linked with the "Port del Comte" alpine ski resort. It is located near the top of the massif (Figure 11), in an area drained by the regional spring M-22, which is the most important resource of the PCM. Given the high karstification degree of the highest parts of the PCM, a hypothetical contamination coming from the ski resort The contamination of groundwater by nitrates in the PCM is mostly related to anthropic activities conducted in aquifer recharge areas. Here, the highest nitrate contents in GW are related to agricultural practices. Other relevant anthropic activities in the study area are restricted to those linked with the "Port del Comte" alpine ski resort. It is located near the top of the massif (Figure 11), in an area drained by the regional spring M-22, which is the most important resource of the PCM. Given the high karstification degree of the highest parts of the PCM, a hypothetical contamination coming from the ski resort would reach the aquifer feeding the spring M-22. Despite this, the impact of the sky resort in M-22 is not relevant at all, at least from the perspective of the NO 3 content in GW. This result stresses the good practices of the ski resort managers in terms of adopting measures to minimize the impact of such activity in the environment.
The dual-isotope diagram δ 15 N NO3 vs. δ 34 S SO4 representing the isotopic composition of atmospheric deposition, soil, fertilizers, and sewage ( Figure A7, Appendix A) for the water samples with both data available (i.e., δ 15 N NO3 and δ 34 S SO4 ; a total of 19 springs) confirms that the main sources of groundwater pollution for springs belonging to cluster A, which are those discharging close to main recharge areas at the top of the PCM, are mainly the atmospheric deposition and fertilizers, and less the mineralization of soil organic matter. Nevertheless, in the case of the regional springs M-22, M-32, and M-43 (Figure 1), the oxidation of sulfides or organogenic sulfur appear to be an additional polluting source as pointed out by the BBM model. Besides, in the case of springs belonging to cluster B, there must be another contribution of sulfate, along with the atmospheric deposition, to explain their isotope content. According to the results of BMM, this source might be the sulfate of geogenic origin, due to the dissolution of Triassic and/or Tertiary gypsum.

Proportional Contribution of NO 3 Sources in GW in PCM
A Bayesian isotope mixing model was prepared using the SIMMR package for R and was run to estimate proportional contributions of the NO 3 source for the 20 spring groundwater samples. The considered sources are the same as the biplots ( Figure 12A) plus atmospheric deposition: (NHF) NO 3 − derived from NH 4 + in chemical fertilizers and precipitation; (NF) NO 3 in chemical fertilizer; (SN) soil organic nitrogen; (S) soil; (M) manure; and (Natm) atmospheric deposition. Figure A6, Appendix A, presents the corresponding outputs separated for each water spring.
The result of the model confirms that, in general, the greatest contribution of nitrates comes from pollution related to anthropic activities carried out in the aquifer recharge areas, mainly the use of fertilizers (NHF), except in a specific case with a notable proportion of manure (M) and sewage (S). Most springs have nitrate concentrations below current standards for drinking water. The only spring that exceeds the reference levels established at 50 ppm is the M-32 spring (with 57.3 ppm). This spring is located downstream of an area of field potato crops. The model indicates that it has an NHF proportion of 32%. The rest of the springs present similar or slightly higher NHF contributions although with lower nitrate concentrations. The second spring with the highest nitrate concentration corresponds to spring M-28 (38.6 ppm). In this case, the origin is clearly influenced by the livestock activity located upstream, presenting a proportion of 39.2% coming from manure and 31.7% coming from sewage.

Conceptual Model for Hydrogeochemical Evolution of GW in the PCM
From the combined analysis of the geological and hydrogeological context and chemical and isotopic data, global hydrogeological and hydrogeochemical conceptual models were interpreted based on the four cross-sections indicated in Figure 1.
The PCM is a high mountain karst aquifer, built upon several thrust sheets of carbonate materials. Precipitation is usually as snow in the highest part of the massif, where bare land abounds, along with the most developed karst forms in the epikarst (Figures 13 and 14). The meteoric water from snowmelt and rainfall infiltrates and recharges the aquifer, mostly as Ca-Cl-HCO 3 type water. Recharged water flows in all directions and discharges through multiple significant springs. The main aquifer of the system is the one associated with the karstified limestones and dolostones of the Tertiary PPEc unit (Figures 1, 13 and 14), which underlies the PEcp1 and PEcp2 units. The four most important springs in the system-named in decreasing discharge rateare M-22 (in Figure 13A), M-25 and M-31 (in Figure 14B), and M-43 (see Figure 1). These springs drain the Tertiary karst aquifer along the syncline axes and were classified into cluster A by Herms et al. (2021) [59]. According to both their recharge elevation zones obtained with the H and O water isotopic content and the 3-D geological structure, these four springs present recharge-discharge pathways, several km long, while presenting a Ca-HCO 3 water type with low TDS (from 122 to 182 mg/L). These characteristics are interpreted as an indicator of a high karstification degree affecting the geological PPEc, PEcp1, and PEcp2 units, which favors large flow rates in both the percolating meteoric water through the unsaturated zone and the GW flow in the saturated zone. The hydrochemical signature of the GW flowing through the karst aquifer is obtained quickly, during the percolation and the first stages of the GW flow phases, as supported by the inverse modeling analysis done with the help of PHREEQC. To illustrate this, Figure 13A shows the conceptual long RDP associated with the regional karst spring M-22. It includes a thick unsaturated zone and a regional water table level at elevations between 1000 and 1100 m a.s.l. Additionally, small and local springs, such as M-08 and M-05 (in Figure 13A), M-15 (in Figure 13B), or M-24 (in Figure 14A), drain the same karst aquifer and have a Ca-HCO 3 hydrochemical water type, with TDS mostly between 99 and 255 mg/L, and SO 4 coming from soil and atmospheric deposition, a composition which is similar to that of the regional springs. Additionally, the local springs may be affected by NH 4 fertilizers and/or manure. The discharge of springs located in the SE part of the study zone, which crosses the limits of the South Pyrenees thrust fault in the front SE of the PCM (right side of the Figure 13B), may be affected by Tertiary gypsum (PExb unit), generated by Ca-SO 4 -HCO 3 to Ca-SO 4 water types.
Below the Tertiary limestone layers, the Garumnian Kgp unit acts as an aquitard, while the Upper Cretaceous Kat, KMca units, the Keuper Tk unit, and the Muschelkalk Tm unit (Triassic) act as local aquifers. These aquifer units drain through small local springs that may have been recharged through the overlying Tertiary carbonate units. The incoming recharge presents an initial Ca-HCO 3 signature, but it changes along the GW flow line by incorporating other solutes from the most soluble evaporite minerals in such local aquifers. There are some springs whose discharge present some hydrochemical special characteristics. Spring M-38 in cluster A ( Figure 14A) interacts with lignite-bearing marls (KMca unit), incorporating sulfate from sulfide minerals (disseminated pyrites) or coal organic sulfur weathering, as is supported by the SO 4 isotope composition in groundwater; spring M-20 in cluster C incorporates Cl by dissolution of halite from the outcropping Keuper materials in the NE part of the massif (Figure 1). Additionally, Ca is dissolved by a reverse ion exchange-like process in weathering of shales, in which Na in dissolution replaces Ca in the terrain matrix. The springs that interact with the Keuper (Tk) unit, which are recharged either in the outcrops of this TK unit or through the geological units overlying it (e.g., Kgp, Kat, KMca), incorporate significant amounts of sulfate by dissolution of gypsum, as happens in springs M-10 (in Figures 13A and 14A) and M-09 ( Figure 14A). Both the S and O from the dissolved sulfate isotope composition and the inverse modeling with PHREEQC support this. Additionally, GW discharge may experience local de-dolomitization, as in the case of M-28 (see Figure 8), which presents the highest content of sulfate (Figure 9), thus inducing the precipitation of calcite, while those with deeper flow lines, such as springs M-30 ( Figure 14B) and M-41 (see Figure 1) of cluster D, incorporate Cl and Na by dissolution of halite as well.

Conclusions
The Port del Comte Massif (PCM) contains one of the most important karst aquifers in the South-Eastern part of the Pyrenees. In this work, hydrochemical and multi-isotope data along with hydrogeological framework information were coupled to characterize the hydrochemical processes driving the hydrogeochemical behavior of this complex hydrogeological system.
In general, the groundwater is dominantly of the calcium bicarbonate and calciummagnesium bicarbonate type, suggesting a dominant calcite dissolution process in agreement with the lithology associated with the Eocene carbonate materials conforming the main aquifer of PCM. The main source of sulfate in GW is the dissolution of geogenic origin from gypsum dissolution from the Eocene-Oligocene Beuda Formation and from Triassic evaporites. Some influence of sulfate from sulfide mineral or coal organic sulfur weathering was also indicated. From the anthropogenic point of view, sulphate from fertilizers seems to play a role in some places around the ski resort. Due to their intermediate values, the role of soil sulphate and sulphate from atmospheric deposition cannot be discarded. Isotope data showed that the source of recharge is precipitation that enters the system along the mountain slopes, favored by the high karstification of the carbonate materials that abundantly crop out in the area. Isotopes also show that dissolved NO 3 in groundwater mainly comes from mineral fertilizers, soil organic nitrogen, and pig manure application to the fields, with at most minor contributions from sewage. As the other high mountain karst systems, the PCM is very vulnerable to pollution. Here, nitrates from agricultural practices represent the main threat to the pristine waters of the aquifer system despite its low significance. Fortunately, the dissolved nitrate concentration in GW is generally low.
The carbonate karstic aquifer of the PCM is a very complex hydrological system developed in a high mountain environment. The multidisciplinary approach allowed the development of a hydrogeological conceptual model of aquifer system functioning, which is coherent with the available information from previous studies, but which is also consistent with the processes driving the hydrogeochemical and isotopic fingerprint of groundwater in this aquifer system. Acknowledgments: We thank the anonymous reviewers for their constructive comments and suggestions which led to a substantial improvement of the paper as well as all the researchers that has collaborated in this research.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Table A1. Summary of water sample types and analysis done in the research project.

Number of Control Points
Total Field Campaigns     Table A4. Chemical characteristics of major constituents for the 10 snow samples (7 natural and 3 artificial) and 2 water ponds. (samples type T-1 correspond to artificial snow (snow gun); T-2 natural snow (inside sky trail); T-3 natural snow (outside sky trail); and T-4 water ponds for artificial snow production.