
Received: 15 July 2016
Accepted: 15 December 2016
Abstract: The Caldes de Boí geothermal waters show important differences in pH (6.5–9.6) and temperature (15.9ºC–52ºC) despite they have a common origin and a very simple circuit at depth (4km below the recharge area level). These differences are the result of secondary processes such as conductive cooling, mixing with colder shallower waters, and input of external CO2, which affect each spring to a different extent in the terminal part of the thermal circuit. In this paper, the secondary processes that control the geochemical evolution of this system have been addressed using a geochemical dataset spanning over 20 years and combining different approaches: classical geochemical calculations and geochemical modelling. Mixing between a cold and a thermal end-member, cooling and CO2 exchange are the processes affecting the spring waters with different intensity over time. These differences in the intensity of the secondary processes could be controlled by the effect of climate and indirectly by the geomorphological and hydrogeological setting of the different springs. Infiltration recharging the shallow aquifer is dominant during the rainy seasons and the extent of the mixing process is greater, at least in some springs. Moreover, significant rainfall can produce a decrease in the ground temperature favouring the conductive cooling. Finally, the geomorphological settings of the springs determine the thickness and the hydraulic properties of the saturated layer below them and, therefore, they affect the extent of the mixing process between the deep thermal waters and the shallower cold waters. The understanding of the compositional changes in the thermal waters and the main factors that could affect them is a key issue to plan the future management of the geothermal resources of the Caldes de Boí system. Here, we propose to use a simple methodology to assess the effect of those factors, which could affect the quality of the thermal waters for balneotherapy at long-term scale. Furthermore, the methodology used in this study can be applied to other geothermal systems.
Keywords: Geothermal system, Secondary processes, Mixing waters, Conductive cooling, Co2, input, Geochemical modelling techniquImportar imagen.
Introduction
Alkaline geothermal groundwaters in granitic rocks have been studied over decades from different points of view as potential geothermal energy resources or as natural analogues of nuclear waste repositories. A number of studies have focused on the hydrogeochemistry and the characterization of equilibrium conditions at depth providing valuable information about these systems (see for example, Michard and Fouillac, 1980; Michard and Roekens, 1983; Fernández and Banda, 1989; Michard, 1990; Alaux-Negrel et al., 1993; Auqué et al., 1996, 1997, 1998; van Middlesworth and Wood, 1998; Druschel and Rosenberg, 2001; Gimeno et al., 2007; Asta et al., 2012).
In this study, the geothermal alkaline system of the Caldes de Boí spa resort has been studied. This spa resort occupies a surface area of 24 hectares and contains a vast extension of gardens where a total of 37 mineral water springs (20 of them have been included in this study) of different chemical compositions can be found. Some of the springs are located in the resort gardens, although the majority of them have been channelled over to the spa’s treatment building. They are thermal waters with emerging temperatures between 15ºC and 52ºC, very low total dissolved solids (0.04–0.35g/L) and pH values 6.5 to 9.6.
The geological setting, the geochemistry and geochemical modelling of the processes occurring in the Caldes de Boí thermal waters (including geothermometrical calculations), have been reported in earlier studies (Corominas, 1978; Albert et al., 1979; Auqué, 1993; Auqué et al., 1996, 1997, 1998; Moreno et al., 1997; Buil et al., 2002, 2006; Asta et al., 2010). According to these studies, the Caldes de Boí thermal waters have the typical features of most thermal alkaline waters present in granitic rocks (Michard and Roekens, 1983; Michard et al., 1986; Sanjuan et al., 1988; Michard, 1990). These thermal waters circulate through a deep hydrological circuit which branches off in different springs in the final stretch. Although the waters show common chemical characters, there are some differences affecting mainly the pH and temperature values. These differences have been attributed to the effect of secondary processes during the thermal water’s ascent (mixing with cold waters and/ or conductive cooling), and to the CO. exchange with soils and atmosphere (Auqué, 1993; Asta et al., 2010 and references therein).
A considerable part of the appeal of many spas rests in the springs’ status as natural phenomena (Atkinson and Davison, 2002) and its maintenance requires a sustainable use of the thermal water preserving the water temperature and the water quality in order to be used for balneotherapy practices. To identify potential threats to the water quality of Caldes de Boí thermal waters over time (e.g. mixing with superficial cold waters) data over more than 20 years have been collected and a classical and simple geochemical methodology similar to that presented by Auqué (1993) and Asta et al. (2010) has been applied. This methodology allows evaluating the effect of the aforementioned secondary processes on the chemistry and temperature of the thermal waters by using different ion-ion plots and geochemical modelling calculations. Furthermore, the geochemical calculations presented here are useful to evaluate the potential precipitation of mineral in the spa installations, which could cause high costs; and to adopt suitable actions to mitigate it. Therefore, the aims of this research are: i) to evaluate the potential effect of additional factors (such as the seasonal climate and the spring geomorphological setting) in the changes observed in the hydrochemistry of the different thermal waters with time; ii) to assess their influence on the thermal waters modifying their composition and temperature, affecting to their quality for balneotherapy and causing mineral scaling; and iii) to propose monitoring practices to maintain the sustainability of the spa resort. With these objectives, reported data of IGME (1984), ENHER (1985), Buil et al. (2002), Asta et al. (2010) together with the data corresponding to a sampling campaign carried out by the authors of this paper in 2008, were studied.
In summary, this research provides a simple methodology to evaluate the vulnerability of thermal waters at a long-term scale when the system is affected by external factors. Furthermore, it could be used as a basis for the evaluation of the exploitation, management and development of thermal water resources.
Geological, Hydrogeological and Climatic Settings
The Caldes de Boí geothermal system is associated to the Hercynian plutonic complex of La Maladeta, which is located in the central sector of the Axial Zone of the Pyrenees (North of Spain; Fig. 1). La Maladeta plutonic complex consists of two main units, the Aneto Unit, in the western part, dominated by sienogranites and the Boí Unit, in the eastern part and dominated by graniodiorites and monzogranites. These two units are separated by a ductil shear zone (Noguera-Ribagorzana milonite zone) and the Caldes de Boí springs emerge through the Boí Unit.
La Maladeta pluton intrudes into Cambrian-Ordovician and Carboniferous sedimentary formations (De Sitter, 1959; De Sitter et al., 1960; Kleinsmiede, 1960; Zandvliet, 1960; Mey, 1965, 1968; Zwart and Roberti, 1976; Vinchon, 1977; Zwart, 1977, 1979). All these sedimentary rocks are metamorphosed by a low-grade epizonal regional

metamorphism (greenschist facies) (Charlet, 1972; Charlet and Dupuis, 1974; Charlet, 1977; Bourke, 1979; Delgado, 1993; Arranz, 1997) and a contact metamorphic aureole, over the regional metamorphism, extends also around the complex (Charlet, 1972, 1977; Charlet and Dupuis, 1974; Bourke, 1979; Delgado, 1993; Arranz, 1997).
The most striking structural feature of La Maladeta massif related to the geothermal springs is its dense fracture network. Apart from the Noguera-Ribagorzana milonite zone that divides the plutonic complex into the two units, there is an important number of fractures and shear zones affecting La Maladeta massif. These fracture systems enhance the subsurface permeability and exert a great influence on the development of this geothermal system. The fractures are mainly related to the Alpine orogeny and show preferred orientations NE-SW and NW-SE (see Fig. 1). This intense network favours a slow deep infiltration of the water through continuous vertical fractures enhancing the development of the geothermal system and acting as the main water flow paths (Fig. 2). The recharge area corresponds to the region of the Gémena lakes, at 2,350m of altitude, and it is located close (at a distance of 3 to 4km) to the emerging zone (Buil et al., 2002). The reservoir is thought to be very deep, around 4,200m, calculated from the temperature at depth obtained by geothermometrical methods and the geothermal gradient of the area (Auqué, 1993; Buil et al., 2002). The residence time of the thermal waters was calculated to be very long (about 16,000 years based on the 14C data; Buil et al., 2002, 2006). The circuit is thought to be U-shaped (Fig. 2) due to the short separation between recharge and discharge areas (3–4km) and to the great depth of the reservoir (4km below the recharge area level). This feature has also been described for other thermal systems in the Pyrenees (e.g.Luchon, Chevalier-Lemire et al., 1990; Cauterets, Soulé, 1990) and in other granitic geothermal systems (e.g. Druschel and Rosenberg, 2001).
The Caldes de Boí geothermal waters undergo conductive cooling during their rise to surface conditions through the granite fractures. This cooling is about 60°C taking into account that the reservoir temperature has been calculated to be around 110°C (Auqué et al., 1998) and the emerging temperature range of the waters is from 15.9ºC to 50.0°C (Fig. I, Electronic Appendix available at www. geologica-acta.com). Apart from the conductive cooling, other processes are believed to have been occurred in the thermal system such as i) mixing with shallow colder water and ii) CO. exchange with soil or atmospheric environments. These processes have been found to affect the pH and the temperature of the thermal waters and, to less extent, their chemical features (Auqué, 1993; Asta et al., 2010).
The study area is geomorphologically located in a narrow glacial valley with steep slopes (Fig. 3). Till deposits cover the bottom and hillsides of the valley. These sediments consist of a basal till material with a Tables I to VI in the Electronic Appendix, available at www.geologica-acta.com). The geochemical data were checked for their appropriate sampling and analytical methodologies. Only samples with in situ determination of pH and temperature and well documented sampling and analytical methodologies, were considered here. Then, the charge balance error for the reported analyses (calculated with PHREEQC; Parkhurst and Appelo, 1999) was also checked and in the majority of the selected spring waters it was lower than ±10%.
The only unpublished data correspond to the more recent sampling campaign, carried out by the authors of this paper, and a summary of the methodology used is described below.

high content of clayey matrix (>50%) and are associated with glaciofluvial deposits. Some periglacial sediments are also present on the slopes. The morphology of the latter favours the presence of landslide deposits such as debris fans and talus slopes that usually cover the glacial deposits. It is remarkable the existence and preservation of fluvial terraces in the valley, mainly constituted by gravels, pebbles and some boulders, often with imbricated disposition, and sandy matrix (Vilaplana, 1983).
The climate of the region is alpine with continental- mediterranean influences. The mean annual temperature is 7.3ºC and the mean annual rainfall is 1,157mm year-1 (Ninot and Ferré, 2008). Most of the precipitation falls in the autumn and spring. The summer is relatively dry (with occasional rainstorms; mean summer rainfall of 345mm year-1; Ninot and Ferré, 2008), as is the winter when snowfalls alternate with long anticyclonic periods (López- Moreno and García-Ruíz, 2004).
Methodology
Sampling and analytical methods
Data compiled by the Spanish Geological Survey (IGME, 1984), by a local hydroelectric company (ENHER, 1985), by Buil et al. (2002), by Asta et al. (2010) and unpublished data from a 2008 geochemical sampling field campaign, were used in this study (see Water pH and temperature were measured in the field and the water samples were taken in HCl-pre-washed polyethylene bottles rinsed, afterwards, three times with double-distilled water. Polyethylene bottles of 1,000ml were used for anions and 100ml bottles for cations analysis. Samples for cations were filtered through 0.1μm (MILLIPORE filters previously cleaned with ultra pure nitric acid) and then preserved with HNO3 to pH less than 2. Water pH was measured using two different portable pH-meters with automatic temperature correction, an ORION 250A and an ORION 8155SC. One calibrated with pH buffers, 4.0 and 7.0, and the other with pH buffers 7.0 and 10.0. Temperature was measured with the probe connected with the pH-meter. The estimated measurement errors for pH and temperature were ≤0.03pH unit and ±0.5ºC, respectively.
Inductively Coupled Plasma Optical Emission Spectrometry (ICP-OES) was used in order to determine the concentrations of major cations (Na, K, Ca) and silica and Inductively Coupled Plasma Mass Spectrometry (ICP-MS) for the determination of Mg concentration. The detection limits were: 0.05mg/L for Na; 0.1mg/L for K; 0.06mg/L for Ca; 0.06mg/L for SiO2 and in the order of 1μg/L for Mg. The error was estimated to be below 5%. Within 24 hours after sampling total alkalinity was analysed by volumetric titration with H2SO4 0.02N and potentiometric determination of final point at pH 4.5. The determination of Cl- and F- concentrations was performed by ion-selective electrodes and SO 2-concentrations were measured with a spectrophotometer using a turbiditic method with Ba. The average measurement error for all the reported analytical results is less than 5%.
Interpretative methods: ion-ion plots and geochemical modelling
Following the methodology used by Auqué (1993) and Asta et al. (2010, 2012) the springs selected to represent non-thermal, thermal and mixed waters were evaluated by the analysis of binary correlation between elements with a conservative behaviour in the studied system (such as Na, SiO2, F, Li). Chloride is usually selected as a conservative mixing tracer but in these waters its concentration is very low and minor analytical or sampling errors may drastically affect calculated mixing trends. Therefore, following the methodology used in studies of this type of alkaline thermal waters (Criaud and Vuataz, 1984; Auqué, 1993; Chae et al., 2006; Asta et al., 2010; Han et al., 2010), sodium has been selected instead as the tracer element for the mixing process. Its concentrations in the thermal and cold end- members are very different, the analytical uncertainty very low, and it behaves as a non-reactive element in the system.

Finally, speciation-solubility calculations and reaction- path simulations have been performed using PHREEQC (Parkhurst and Appelo, 1999) with the WATEQ4F thermodynamic database (Ball and Nordstrom, 2001). The solubility data for kaolinite proposed by Michard (1983) have been included in this database, as it has been successfully used in the geochemical modelling of thermal waters in many different crystalline systems (e.g. Michard et al., 1986; Auqué et al., 1998; Buil et al., 2006; Gimeno et al., 2007; Sanjuan et al., 2010; Asta et al., 2012). These calculations were aimed to identify and quantify the different secondary processes that affect the system and they have been done in closed system conditions (e.g. with respect to CO2) except when indicated otherwise.
Climatic information and water balance calculation
The meteorological information was obtained from the meteorological stations at Caldes de Boí and Pont de Suert (located at 50km of the studied area). A time period of 20 years (1983–2003) was selected and the monthly precipitation and temperature data were collected from the Meteorological Service of Catalonia data base. Since there were not available precipitation data for 1997 in Caldes de Boí, the data from Pont de Suert were used instead after checking that the differences in the monthly rainfall averages for both locations during a period of 58 years were less than 16%.
The Thornthwaite-Matter method (Dunne, 1978) and the evapotranspiration Thornthwaite formula (Thornthwaite, 1948) were used to calculate the water balance. The input variables are: mean monthly temperature (in degrees Celsius), monthly Total Precipitation (TP, in millimetres), Soil Moisture Capacity (SMC) and latitude (in decimal degrees) of the location of interest.
Most of the recharge area corresponds to granite outcrops with very little soil development. When soils are present they correspond to leptosols (Jones et al., 2005), shallow and sandy soils with low water retention capacity and high permeability (high infiltration potential). Since no prior studies on soil morphology exist for this study area, the soil-moisture storage capacity was estimated using general information about this type of soil (Lal and Shukla, 2004; Gupta et al., 2006). For these soils we estimate that soil moisture capacity is between 20 and 50mm. The output of the model consists of monthly water balance including monthly potential and Actual EvapoTranspiration (AET), soil moisture storage and surplus. Actual evapotranspiration is derived from potential evapotranspiration, total precipitation and soil moisture storage. In this model, the soil is assumed to have a maximum soil moisture capacity. Then, moisture is added to or subtracted from the soil, depending on whether the monthly precipitation is greater or less than the monthly potential evapotranspiration. If TP exceeds Potential Evapotranspiration (PET), then AET is equal to PET and the water in excess of PET replenishes the SMC and the excess of precipitation is assumed to contribute to the water surplus. If the sum of TP and SMC is less than the PET, then a water deficit is calculated as PET-AET. Studying the evolution of the water surplus or deficits, it is possible to estimate recharge during dry or humid periods. Infiltration is used to estimate recharge taking into account the climatic, topographic and soil characteristics.
Results and discussion
Identification and verification of the secondary processes
Figures 4 and 5 present several ion-ion plots in which, as in earlier studies (cf Criaud and Vuataz, 1984; Auqué, 1993; Chae et al., 2006; Asta et al., 2010, 2012; Han et al., 2010), sodium has been selected as the main conservative geothermal tracer in this system. A very good correlation is shown between conservative geothermal tracers for this system including SiO2, F, Li and Na. This suggests a mixing process between two end-members with different chemical characteristics and corresponding to a thermal end-member (represented by Estufa, Termas and Baños springs) and a cold end-member (represented by Ferro spring). This mixing relationship appears to be time invariant. The mixing fraction for each spring can be calculated by its position on the mixing line (see Tables II to V in the Electronic Appendix). Although this mixing process occurs throughout the time, the mixing proportion can be different at different springs. This will be analysed below.
The presence of secondary processes has also been checked by analysing the relationship between the concentration of sodium (as the conservative element in this case) and other key physico-chemical parameters such as temperature and pH. Figure 6 shows the evolution of sodium concentration with the cooling intensity. Formally, the cooling intensity can be expressed in terms of enthalpy following the approach by Fournier (1979). However, in the temperature range of interest (below 50ºC), there are not meaningful differences between temperature and enthalpy and, thus, the cooling intensity has been expressed as the difference between the maximum temperature in the thermal end-member, T0, and the emergence spring temperature, T, in each sampling campaign (from zero for the thermal end-member (Estufa spring) to almost 40ºC for the cold end-member (Ferro spring)).
If no additional secondary cooling processes were acting in the system, sodium content would show a linear and negative correlation with the cooling intensity due to mixing. Indeed, as observed in Figure 6, there are some waters placed on the mixing line that are clearly affected by mixing with cold shallow water. In contrast, the springs situated to the right of the mixing line (e.g. Santa Lucía, Tartera L, Tartera R, Salenca, Pompeio, Termas and Bosque) are affected by additional cooling which has been interpreted as the result of conductive cooling along longer water flow paths. When the location of these points is extended in the horizontal line to intercept the mixing line the degree of additional cooling can be estimated.
Three trends can be observed by analysing the evolution of pH with respect to the cooling (Fig. 7) from all sampling campaigns: i) a slight and almost linear pH increase as conductive cooling increases; ii) a more noticeable pH decrease as cooling increases, outlined by the springs only affected by mixing (compare with plots shown in Figure 5) and iii) a decrease of pH from the hot end-member to the pH values measured in one of the warmest waters (Termas).
With respect to this last trend, and as stated above, Estufa and Termas have both high temperatures (48.8ºC to 52ºC and 45.8ºC to 50.1ºC, respectively), low conductive cooling, and almost identical chemical compositions, but different pH values (8.8 to 9.34 and 7.8 to 7.99, respectively; Fig. 7). This behaviour has been interpreted as the presence of an additional process affecting only the pH of the solutions without affecting




the rest of the chemical characteristics of the fluid. Carbon speciation-solubility results (Fig. 8) indicate that Estufa (thermal end-member) shows a noticeable disequilibrium with respect to the atmospheric CO2 partial pressure (log pCO2 between -4.01 and -4.92) while the pCO2 obtained for Termas is very close to equilibrium (log pCO2 between -3.03 and 3.21). This fact suggests that, although Termas could have had an original CO2 partial pressure similar to Estufa, when it comes in contact with the atmosphere or soil environment it is equilibrated with its CO2 content (incorporating CO2 into solution). This does not happen to Estufa before the emergence point. This CO2 input could be the responsible for the decrease in pH compared with the Estufa pH value. Below we evaluate all of these hypotheses using geochemical modelling.
Geochemical modelling
Following the procedure described in Asta et al. (2010), simulation of the different secondary processes affecting the system has been performed with PHREEQC (Parkhurst and Appelo, 1999). The simulations include the pH evolution along three processes: mixing, cooling and CO2 exchange and finally their superposition. The data used as end-members correspond to each sampling campaign, except for the case of the 1997 sampling, when the cold end-member used corresponds to Ferro data of 2003 or in the case of 2008 campaign when the hot end members, Estufa and Termas, were not sampled and the data from the 2003 campaign have been used instead.
Figure 9 shows the evolution of pH with respect to cooling (T0-T). The calculated evolution of pH during mixing between the thermal (some of them affected also by CO2 input) and cold end members is shown as solid and discontinuous lines connecting the different end members used for each mixing process. The pH values are very similar (or slightly higher) for a wide range of mixing proportions of the hot end-member. Only when the cold end-member is in a proportion greater than 40% (0.6/0.4 mixing proportion) the pH decreases towards the value observed in the cold water. This is the result of the important pH buffering capacity of the thermal water due to its high pH values and high concentrations of CO 2- and H3SiO - and the presence of unusually active acid-base pairs (HCO -/CO 2- and mainly H SiO /H SiO -) in these alkaline waters acting as homogeneous buffers (Asta et al., 2010). In contrast to the Estufa-Ferro mixing simulation, the result obtained using the less alkaline hot end-member, Termas, does not show the buffered region in any of the sample sets, as its CO 2- and H SiO - contents are much lower. The comparison of the pH and cooling values of the sampled waters with these paths (Fig. 9) shows that some springs fit well in the theoretical pH lines due to mixing between the end members (Estufa and Ferro or Termas and Ferro) (see for example Baños, Santa Lucía, Tartera L, Canem 2 in 1997 data, or Santa Lucía, Titus 1, Titus 2 and Avellaner in 2003 data; Fig. 9). However, neither these spring waters consistently fit in the mixing lines, nor they are mixed in the same mixing proportion. The rest of the springs, which do not follow the end-member mixing lines, show pH values that must have been affected, to some extent, by other geochemical processes.
Other acid-base pairs, such as, boron, ammonia and sulfide (H3BO3/H2BO3 - , NH4 +/NH3, and H2S/HS- ) could be active in thermal alkaline waters buffering the pH. However, those species are in very low concentrations in the Caldes de Boí thermal waters. For example, boron concentrations in the 1983 sampling (IGME, 1984) were between 0.016 and 0.086mmol/L and in the 1985 sampling (ENHER, 1985) between 0.018 and 0.042mmol/L with the highest values found in the Estufa spring. Using these values in the speciation-solubility calculations (data not shown) the results showed that boron and ammonia concentrations are two orders of magnitude lower than those of the H4SiO4 and H3SiO4 - species and almost one order of magnitude lower than the concentration of the OH- species. In addition, to check the effect of the presence of these species on the pH, mixing calculations including them were performed with the 1985 data. The results showed that the buffering contribution of those acid-base pairs would be of minor importance and were not considered in this study.
No data were available on the total dissolved sulphide species for the Caldes de Boí thermal system.

The presence of some sulphide species is possible in reduced thermal waters (e.g. deduced from the available Eh values), however, the good charge balances obtained for the Estufa spring (better than ±3%) suggest that the concentration of total dissolved S-2 is low. With respect to the cold end-member, its redox character (available Eh values around +200mV) indicates the absence of this component. On another way of reasoning H2S/HS- is an acid-base pair and also a redox pair, and its behaviour would be mainly controlled by redox-related parameters during mixing processes. For example, Criaud and Vuataz (1984) found that, during the mixing processes at the Luchon thermal system (an alkaline thermal system in the French Pyrenees), dissolved sulphide in the thermal end- member was oxidised to sulphate by the oxygen dissolved in the cold shallow end-member. Therefore, the possible buffering capacity H2S/HS- would be very limited.
Furthermore, to evaluate the impact of the uncertainty in the Ferro spring composition (with charge balance errors up to -12.37%) a sensitivity analysis has been performed in the mixing calculations using the 2003 and 2005 data sets (with the highest charge balance errors). The mixing simulations have been made: i) using Ferroas the cold end-member, with their actual compositions, and ii) forcing the waters to have a balance charge of 0% changing the Mg concentration. The results showed no significant differences supporting the mixing model presented in this study.
The second PHREEQC simulation was designed to assess the effect of conductive cooling by means of a progressive decrease in temperature over the thermal solutions Estufa and Termas in successive steps of 5°C of decrease. Additional simulations were also considered to show the effect of the conductive cooling of waters already mixed in different proportions (0.7/0.3 and 0.5/0.5). The effect on pH is a linear and progressive increase in pH as the temperature decreases (straight lines in Figure 9). The results show a pH increase rate of 0.012–0.013pH units per degree centigrade (of temperature decrease), which is in good agreement with the pH increase rate (0.013 units per degree centigrade) reported by previous authors (Boulègue, 1979; Michard and Fouillac, 1980; Auqué, 1993).
The effect of conductive cooling and mixing between hot and cold solutions on pH is clearly different when the cold end-member is in a proportion greater than 40%. As

already reported by Asta et al. (2010), when the hot water dominates the mixing proportion, the conductive cooling effect on pH may be undistinguishable from the cooling effects due to mixing (see for example in Figure 9, Tartera L in the 1997 and 2003 samplings). This is because the pH variation is very small at the beginning of both processes. When the proportion of cold water increases in the mixing, the pH decrease dominates over the opposite trend of pH increase due to conductive cooling. This effect is also seen when the hot end-member is a less alkaline water (Termas spring), and the buffering effect of pH during mixing is much lower. In that case the effect of decreasing the pH due to mixing is higher.
The fact that the pH values do not fit exactly on the theoretical mixing line indicates that still some additional process must be considered. As described above, CO2 partial pressure (pCO2) calculated for these waters range from logpCO2= -5.2 (Salenca in the sampling campaign of 1983) to logpCO2= -1.90 (Titus 1, in 2008) (see Figure 8 and Tables II to V), and it is expected that waters with low CO2 partial pressure will tend to equilibrate with the atmosphere or the soil system (to incorporate CO2 and therefore, to decrease pH) when they reach surface. This CO2 input can occur in the hot water or after its mixing with some proportion of cold waters. The clearest example shown in Asta et al. (2010), and confirmed for the rest of the sampling campaigns, corresponds to Estufa. The pCO2 value in the thermal end-member Estufa is in clear disequilibrium with the atmosphere, while is close to equilibrium in Termas (Fig. 8). This suggests that Termas, although originally could have had a value similar to Estufa, has experienced the addition of atmospheric/soil CO2. This CO2 input to the Termas thermal water would be responsible for its pH decrease, in contrast to the constant pH of Estufa, which is unaffected by additional CO2 input. Other possible processes affecting the pH value, such as calcite dissolution/precipitation reactions, are not probable (see further details in Asta et al., 2010).
Additional PHREEQC simulations were performed to assess the effect of the CO2 input on the thermal water geochemical composition. To test the ability of the end members to be the sources of mixing, the hot end-member from each sampling campaign was subjected to additional CO2 to generate initial solutions with different pH values. Then, a mixing process with the cold end-member in different proportions is imposed on each initial solution. The final pH values obtained in the simulations are shown in Figure 9 as dotted lines. Higher values of CO2 partial pressures result in lower pH values and then, in lower pH buffering capacity of the waters.
The comparison between the actual values and the results of this simulation of CO2 input confirms the hypothesis proposed above for the Termas spring but also for other waters. This is the case of Canem 2 and Bou in the 1997 sampling and Titus 1 and 2, and Avellaner in the 2003 data (Fig. 9) where the cooling is only due to mixing and their pH values are controlled by mixing with the thermal water Termas (whose CO2 has already been modified). Another possibility is that mixing with the thermal water (Estufa), took place first, and then the mixed waters were affected by the CO2 input; see for example Canem 1, Bou and Avellaner in the 1983 sampling or Tartera L and R, Salenca, Pentanca R and L, Pompeio, Canem 1, Titus 1 and Titus 2 in 1985, or Canem 1 and Besort in 2003, and Titus 1 in 2008.
These results indicate that the combination of mixing, conductive cooling and CO2 input justifies the chemical composition of these waters and the fact that having a similar chemistry, their pH values are quite different. Some springs fit perfectly well with a single mixing process, either between the original end members, Estufa and Ferro, or between a modified hot end-member (Termas) and the same cold end-member (Ferro). Conductive cooling produces a pH increase that is undistinguishable from the cooling due to mixing when the hot end-member dominates the mixing proportion. The input of CO2, finally, has been shown to be present in the system producing a pH decrease, which is evident in the springs located below the theoretical mixing lines.
Temporal evolution of the secondary processes
According to the results shown in the previous sections, the mixing between a cold and a thermal end-member is a generalised process that affects most of the spring water system over time (Fig. 10). Mixing is apparent from the good correlation between conservative elements (Figs. 4; 5). The geochemical modelling performed has confirmed that the Caldes de Boí system is also affected by other processes including: conductive cooling and CO2 exchange. These processes have been found to vary over time affecting the spring waters with different intensity. Figure 10 is a summary sketch of the processes and their quantification in the sampling campaigns considered in this study.
The hot and cold end members (Estufa and Ferro, respectively) are relatively constant in temperature and pH over time and have been considered as not affected by the secondary processes. The rest of the spring waters, however, are affected by mixing, conductive cooling and CO2 input to different extent (Fig. 10). The study of these processes over time has shown some interesting behaviours that have helped to understand the external factors on the physico-chemical characters of this geothermal system.

There is a group of springs (Canem 1, Bou, Avellaner 1, Titus 1, Titus 2 and Besort) that are affected mainly by mixing and CO2 input in all sampling campaigns, causing slight variations in temperature and pH. Only in the sampling of 1985 they are also affected by conductive cooling. These variations are due to the fact that the intensity of the studied secondary processes could vary seasonally.
Salenca spring is affected by mixing with similar proportions of the thermal and cold end members (from 0.7/0.3 to 0.8/0.2) and the conductive cooling affects the spring waters in a similar degree in all campaigns (from 12.3ºC to 17ºC of conductive cooling). Only in the 1985 and 2003 samplings Salenca spring water is also influenced by CO2 addition, which could explain the lowest pH values found in those samples. This input of CO2 to the solution could come from the soil or biogenic processes, however, some CO2 addition from the atmosphere during the sampling cannot be ruled out.
From 1983 to 1997 Pompeio spring is affected by mixing with very different proportions of the thermal and cold end members (from 0.3/0.7 to 0.9/0.1). Additionally Pompeio experiences conductive cooling in 1985 and 1997, and CO2 input in 1983 and 1997. In these two samplings Pompeio shows similar pH and spring temperatures.
Other springs such as Tartera R and Tartera L show also high temporal variability with respect to their pH and temperature. This variability is determined by the different processes and by their different intensity. For example, Tartera R shows features very similar to the thermal end- member in 1983, though with lower temperature, suggesting the effect of conductive cooling. In 1985, this spring is influenced by mixing and CO2 addition. And, finally, in 1997 and 2003 the spring shows the effect of all these processes. In the case of Tartera L, it experienced mixing, conductive cooling and CO2 addition in 1983 and 1985, whereas in 1997 the conductive cooling is the only process that plays and important role (displaying similar composition to the thermal end-member). Finally, in 2003, conductive cooling and mixing take place as the main processes.
The high variability in the intensity of the secondary processes is reflected in the physicochemical characteristics of the thermal waters; see for example the evolution of the temperature and pH of Tartera R and Tartera L through time (Fig. 10). In contrast, other spring waters, such as Santa Lucía, Baños, Bosque or even Termas, have relatively constant temperatures and geochemical characteristics, suggesting that they are not affected by temporal variations in fluid mixing or CO2 inputs.
Seasonal climate effect on the water chemistry
Temporal variation in the recharge of the shallow aquifer could affect the extent of the mixing, as well as the rates of conductive cooling of the spring water. To evaluate the potential influence of the climate on the geothermal system, the water balance of the Caldes de Boí area was carried out. Not all the data necessary to perform a complete water balance (e.g. daily climate data, soil analysis, etc.) were available, but this was sufficient to assess if the samplings were carried out in dry or rainy seasons.
The results shown in Figure 11 indicate that the 1983 and 1997 sampling campaigns were carried out during relatively arid conditions without soil water surplus. In contrast, the 1985 sampling clearly corresponds to a period of elevated precipitation with water surplus and with a remarkable water surplus in the soil in the months previous to this sampling.
The extent of mixing appears to be affected by temporal variations in the recharge of the shallow aquifer. There is a group of springs (e.g. Tartera R, Santa Lucía, Pompeio) that are mainly affected by conductive cooling in the dry periods whereas in the rainy seasons the same springs are also affected, although in a different extent, by mixing processes with cold waters of the system. This is due to the fact that during the rainy season is when the shallow aquifer is recharged. Significant rainfall influences also the ground temperature (van Manen and Wallin, 2012), which favours the conductive cooling process. This effect is observed in the waters sampled in 1985 (see Fig. 10) when, conductive cooling affects some springs that are not affected by it in the rest of the samplings (e.g. Titus 1&2, Avellaner).
Large differences between the pH trends amongst different sampling campaigns are observed, with the extreme case of 1985 data when the pH trend is not correlated with any other parameter due to the instability caused by the rainy period (Fig. 9). In contrast, in the samplings carried out during dry periods (1983 and 1997), the three trends indicated in Figure 7 are clearly observed.

Comparing the values of pH with temperature, the variations observed in the secondary processes over time could be due to the hydrological impact of variations in recharge on the geothermal system, mainly through its direct effect on the physico-chemical characters of the cold end-member, representative of the shallow aquifer. At the same time, the shallow aquifer has influence on some factors related to the effects of the thermal energy attenuation in the system (e.g. higher contribution of cold waters or oscillations in the water table of the shallow aquifers).
Dry seasons correspond to periods in which water deficit exists and a progressive discharge of the shallow aquifer takes place. In contrast, in periods of high recharge there is a predominance of infiltration and therefore a recharge of the shallow aquifer. These high infiltration events increase, in general, the extent of the mixing process. The pH trends are less obvious than those observed during dry periods.
Influence of the geomorphological setting
Despite the fact that climate exerts an important effect over the physico-chemical features of the geothermal waters, rainy-dry periods do not affect all spring waters to the same extent. The differences could be explained by the distinct geomorphological and shallow hydrogeological setting of the springs (Figs. 2; 3).
The springs are located over four different types of superficial deposits: i) scree deposits (Tartera R and Tartera L), ii) slope-attached glacial sediments or till (Estufa, Termas and Santa Lucía), iii) fluvial terraces (Canem, Bou, Avellaner and Titus), and iv) debris fan deposits (Salenca) (Fig. 3). In high mountain areas, the glacial and fluvial sediments accumulated at the bottom of the valley may host a shallow alluvial aquifer and only a minor subsurface flow may be present in the weathered surficial bedrock layer and in the highly permeable slope and till deposits attached to the hillside. Thus, it is expected that the deep thermal waters discharging here are covered by fluvial terraces and pass through the alluvial aquifer of the Noguera de Tor River, increasing the mixing between the thermal and the cold waters. This is the case of Canem, Bou, Avellaner and Titus springs. On the other hand, Estufa, Termas, Santa Lucía and Tartera R and Tartera L are found around the limit between the slope-attached till and scree deposits, and therefore, in this setting the deep thermal waters will be only mixed with the hillside subsurface flow waters limiting the possibility to affect significantly the composition of the deep thermal waters.
Salenca, in the lowest part of a debris fan close to the riverbed, shows intermediate features that do not change over time. This suggests that Salenca is not influenced by the shallow aquifer hosted in the alluvial deposits of the Noguera de Tor River and its geochemistry may be influenced by a slightly different flow path followed by the deep geothermal waters.
Practical applications of the mixing model
Our results showed that the geomorphological setting of the springs and the temporal variation in the recharge of the shallow aquifer affect the extent of the mixing, which could cause a decrease in water temperature and a dilution of the thermal water.
Balneotherapy is defined as the use of baths (in tubs or pools) containing thermal and/or mineral water from natural springs or drilled wells with a natural temperature of at least 20ºC and/or mineral water with a total mineral content of at least 1g/L (Gutenbrunner et al., 2010). Then, the increase of cold water fractions in the thermal water would affect the temperature and also the composition and quality of the water. A method to evaluate the impact on the recharge in the aquifer is using mixing models. These models have been used in previous studies to quantify exchanges intra- and inter-aquifer (Gómez et al., 2014), to estimate the degree of interaction between surface and groundwater (Oyarzún et al., 2014), and the relative inflow contributions to the alluvial aquifer (Martinez et al., 2017).
The impact of the cold end-member temperature during the mixing process was tested decreasing the temperature of the cold end-member from 12.3ºC to 1ºC similar to summer and winter temperatures, respectively. The results show no dramatic differences between the mixing using this set of temperatures (see Fig. 12). In general, mixing cold water fractions of 0.6 or more (corresponding to thermal/non-thermal water ratios of 0-0.4) would decrease the thermal water T to values lower than 20ºC. In addition to that, the introduction and mixing of oxygenated cold water could quickly oxidize the sulphide present in the water and dilute the concentrations of other components that provide the water their balneotherapeutic properties (e.g. silica, sodium, chloride). This indicates that thermal/ non-thermal mixing ratios lower than 0.4 would threaten the thermal water quality for the spa purposes.
Differences in the CO2 inputs (or loss), and, therefore, in the pH values have been observed over time in the different springs due to external factors. Although the variations observed are not sufficiently important to affect the water used for baths they could have a critical impact on the spa installations, causing the formation of mineral deposits (e.g. calcite or kaolinite). Mineral scaling in wells

and other installations is a very common process associated with thermal waters (Kele et al., 2015; Todorovic et al., 2016; Asta et al., 2017) and could cause important costs to the spas. The effect of CO2 concentrations in the cold end- member, more susceptible to changes due to exchange with the atmosphere or CO2 from the soil, during the mixing was tested by increasing and decreasing the pCO2 from
In contrast to the calcite results, kaolinite saturation indices are affected by increasing and decreasing the pCO2 and their consequent pH variations (Fig. 13B). The presence of kaolinite precipitates have been detected in springs of other alkaline thermal waters in the Pyrenees (e.g. Boulègue et al., 1981; Couturier et al., 1984) and the obtained results would suggest the possible precipitation of this mineral phase during mixing processes with a high pCO2 cold end-member (Fig. 13B). The importance and implications of this precipitation process would need further studies.
An important component in the management of the thermal water resources is the understanding of the geochemical processes affecting the water and their effects on its properties. In this study we have assessed temporal and spatial variations in spring chemistry as the result of natural external processes. This allowed us to identify the processes that could threaten the quality of the thermal waters. Furthermore, natural and anthropogenic variations associated with groundwater extraction from spas or changes in irrigation practices could increase their mixing with non-thermal water.
The results show that the quality of the thermal water is greatly influenced by the introduction of shallow cold water. Therefore, it is important monitoring changes in water composition (e.g. decrease of Na, F, Li, which behave conservative in the Caldes de Boí system) and temperature over time. This could be used as an early warning for potential problems related to an increase in the water recharge that will cause the decrease in the temperature and reduce the quality of the thermal waters.
Our mixing model showed that cold-water proportions over 60% could be critical for the water properties. Furthermore, the mixing with cold waters could modify the quality of the thermal waters if superficial waters are affected by contamination.
Conclucions
Classical geochemical calculations and geochemical modelling have been used to study the hydrogeochemical evolution of the Caldes de Boí geothermal system over more than 20 years. The results suggest that the hydrogeochemistry of the system could have been mainly determined by the superimposed effects of the following irreversible processes: i) mixing between a cold and a thermal end- member; ii) conductive cooling and iii) CO. input from external sources (atmosphere and soil environments). These processes affected the spring waters with different intensity in the different springs and over time. This suggests that this geothermal system is a dynamic system affected by additional external factors that cause variations in the intensity of the secondary processes and, consequently, in the physicochemical features of the spring waters.
The observed differences in the intensity of the secondary processes over time may be attributed, to some extent, to the effect of climate and indirectly to the geomorphological/hydrogeological setting of the different springs. The mixing process is, in general, greater in rainy seasons, probably due to the increased infiltration and recharge of the shallow aquifer. At the same time, the rainfall exerts influence in other factors related to the thermal energy in the system, such as the decrease in ground temperature promoting the conductive cooling process. The geomorphological/hydrogeological setting of the springs determines the thickness and features of the shallow saturated zone and, therefore, the possibility and the intensity of the mixing between the surficial cold waters and the deep thermal water. To explain the observed variations over time in the studied geothermal waters future research should focus on the study of possible additional external effects in collaboration with the spa, on the detailed characterisation of the thickness of surficial deposits and their saturated zone or on the presence of fractures using geophysical methods. These parallel studies could provide enough information to confirm the hypotheses derived from the geochemical analysis

The results of this study show that the quality of the thermal water is significantly influenced by the introduction of shallow cold water. The model shows that cold-water proportions over 60% could be critical for the water properties of the studied system to be used for balneotherapy practices causing a decrease in temperature (lower than 20°C) and in the concentration of important elements that provide the waters their balneotherapeutic properties. In addition, mixing with cold waters would increase the vulnerability of the thermal waters to be affected by superficial contamination. Finally, it is worth mentioning that the results of this study are important because understanding the changes in the groundwater composition due to mixing, temperature variations and/ or CO2 exchange could be used for the implementation of preventive actions in the Caldes de Boí geothermal system. Regular chemical monitoring of the thermal water and temperature should be done in order to establish natural fluctuations and any longer term trends, which may indicate increase in the introduction of cold water or anthropogenic impacts on the thermal water.
The use of the simple methodology presented in this work allows assessing the effect of external factors, which could affect the quality of the thermal waters for balneotherapy, at long-term scale. This methodology has been tested for a long period of time and can be applied to the study of other low enthalpy geothermal systems, allowing an accurate prediction of their future behaviour and a realistic planning of their future management. For example, spring compositional changes could indicate an increase of the cold meteoric recharge suggesting an overexploitation of the geothermal resources, which could be critical for bathhouses and spas, where specific hot springs temperatures are required. Furthermore, the methodology used in this study can be used to provide some insight for the use of these types of geothermal systems as natural analogues for the emplacement of nuclear wastes. For example, to evaluate the long term evolution of the artificial hydrotermal systems developed during the emplacement of spent fuel (e.g. Druschel and Rosenberg, 2001), the evolution of the alkaline plumes promoted by the interaction of groundwater with cementitious materials (Savage, 2011) or to simulate the effect of mixing processes, especially important in some storage concepts (e.g. the Finnish or Swedish cases; see Gimeno et al., 2014).
Acknowledgements
The economical support of the Aragon Government (DGA) through their program for financing research in Consolidated Groups is acknowledged. We wish to express our gratitude to the staff of the Caldes de Boí resort for their help during the sampling campaigns. We are also very thankful to Enrique Oliver, from the Earth Sciences Department of Zaragoza University, for his technical assistance. Dr. Asta and Dr. Galve have received economical support from the Spanish Ministry of Science and Innovation through a Research Contract from the “Juan de la Cierva Subprogram”. Dr. Acero has received the economical support from the European Commission through the Marie Curie Research Fellowship Programme. The constructive comments of an anonymous reviewer and the assistance of the Editor, Dr. Jordi Cama, have considerably improved the original manuscript and are gratefully acknowledged.
References
Albert, J.F., Corominas, J., París, C., 1979. El estudio hidrogeológico de los manantiales y su aplicación geológica: caso de las aguas termales, carbónicas y sulfhídricas de Cataluña. Acta Geologica Hispánica, 14, 391-394.
Alaux-Negrel, G., Beaucaire, C., Michard, G., Toulhoat, P., Ouzounian, G., 1993. Trace-metal behaviour in natural granitic waters. Journal of Contaminant Hydrology, 13, 309-325.
Arranz, E., 1997. Petrología del macizo granítico de La Maladeta (Huesca-Lérida): estructura, mineralogía, geoquímica y petrogénesis. Doctoral Thesis. Universidad de Zaragoza, 319pp.
Asta, M.P., Gimeno, M.J., Auqué, L.F., Gómez, J., Acero, P., Lapuente, P., 2010. Secondary processes determining the pH of alkaline waters in crystalline rock systems. Chemical Geology, 276, 41-52.
Asta, M.P., Gimeno, M.J., Auqué, L.F., Gómez, J., Acero, P., Lapuente, P., 2012. Hydrochemistry and geothermometrical modeling of Panticosa geothermal system (Spain). Journal of Volcanology and Geothermal Research, 235-236, 84-95.
Asta, M.P., Auqué, L.F., Sanz, F.J., Gimeno, M.J., Acero, P., Blasco, M., García-Alix, A., Gómez, J., Delgado-Huertas, A., Mandado, J., 2017. Travertines associated with the Alhama- Jaraba thermal waters (NE, Spain): Genesis and geochemistry. Sedimentary Geology, 347, 100-116.
Atkinson, T.C., Davison, R.M., 2002. Is the water still hot? Sustainability and the thermal springs at Bath, England. London, Geological Society, 193 (Special Publications), 15-40.
Auqué, L.F., 1993. Estudio de los sistemas geotermales en Aragón. Pautas de especiación y reacción aplicadas a la modelización de sistemas de baja-media entalpía. Doctoral Thesis. Universidad de Zaragoza, 509pp.
Auqué, L.F., Mandado, J., Gimeno, M.J., López, P.L., Gómez, J., 1996. Los sistemas geotermales del Pirineo Central. I. Caracteres geoquímicos y fisicoquímicos de los manantiales termales. Estudios Geológicos, 52, 161-173.
Auqué, L.F., Mandado, J., López, P.L., Gimeno, M.J., 1997. Los sistemas geotermales del Pirineo Central. II. Resultados de la aplicación de técnicas geotermométricas. Estudios Geológicos, 53, 45-54.
Auqué, L.F., Mandado, J., López, P.L., Lapuente, P., Gimeno, M.J., 1998. Los sistemas geotermales del Pirineo Central. III. Evaluación de las condiciones en profundidad y evolución de las soluciones hidrotermales durante su ascenso. Estudios Geológicos, 54, 25-37.
Ball, J.W., Nordstrom, D., 2001. User’s manual for WATEQ4F with revised thermodynamic database and test cases for calculating speciation of major, trace and redox elements in natural waters. U.S. Geological Survey, Water-Resources Investigations Report, 91-183.
Boulègue, J., 1979. Formation des eaux termales sulfureés des Pyrénées Orientales. Origine du Soufre. Géochimie du fer et du cuivre. Journal Francais d’Hydrologie, 10, 91-102.
Boulègue, J., Fouillac, C., Michard, G., 1981. Dêpots minéraux à l’émergence des sources thermales sulfurées sodiques des Pyrénées Orientales. Bulletin de Mineralogie, 104, 661-668.
Bourke, D., 1979. Etude Géologique de la terminaison orientale du Massif de la Maladeta et de sesabords, región d’Espot (province de Lérida, Pyrénées espagnoles). Doctoral Thesis. Lille University, 69pp.
Buil, B., García, S., Lago, M., Arranz, E., Auqué, L., 2002. Estudio geoquímico de los procesos de interacción agua-roca sobre sistemas geotermales de aguas alcalinas en granitoides. Madrid, ENRESA, Publicación Técnica 02, 246pp.
Buil, B., Gómez, P., Turero, M.J., Garralón, A., Lago, M., Arranz, E., De la Cruz, B., 2006. Factors that control the geochemical evolution of hydrotermal systems of alkaline water in granites in Central Pyrenees (Spain). Journal of Iberian Geology, 32, 283-302.
Chae, G.T., Yun, S.T., Kim, K., Mayer, B., 2006. Hydrogeochemistry of sodium-bicarbonate type bedrock groundwater in the Pocheon spa area, South Korea: water– rock interaction and hydrologic mixing. Journal of Hydrology, 321, 326-343.
Charlet, J.M., 1972. Étude géologique et petrographique du Massif Granitique de La Maladeta (Pyrénées Centrales Espagnoles Doctoral Thesis, Polytechnic University of Mons, 115pp.
Charlet, J.M., 1977. Le métamorphisme au contact des granitoïdes entre les vallés de l’Esera et de la Noguera Ribagorzana (Pyrénées Centrales Espagnoles). Annales de la Société Géologique du Nord, 97, 165-177.
Charlet, J.M., Dupuis, C., 1974. Observations nouvelles dans la Massif de la Maladeta. In Actas del VII Congreso Internacional de Estudios Pirenaicos (La Seo de Urgel, 1974): 37-38.
Chevalier-Lemire, G., Pigassou, R., Rigaill, R., Vilmues, T., 1990. Étude des variations naturelles du debit des sources thermales à Luchon (Haute-Garonne, France) par modèle hydrologique global pluies-débits. Hydrogéologie, 4, 287-296.
Corominas, J., 1978. Condiciones hidrogeológicas de los manantiales sulfhídricos de Cataluña. Acta Geologica Hispánica, 13(1), 26-30.
Couturier, Y., Michard, G., Sarazin, G., 1984. Constantes de formation des complexes hydroxydés de l’aluminium en solution aqueuse de 20 ºC a 70 ºC. Geochimica et Cosmochimica Acta, 48, 649-659.
Criaud, A., Vuataz, D., 1984. Etude géochimique et géothermique des eaux sulfurées sodiques de Luchon, Pyrénées. Rapport du Bureau de Recherches Geologiques et Minières (BRGM) Rapport 84 SGN 384 IRG, 61pp.
De Sitter, L.U., 1959. Geological map of the Central Pyrenees, 1/50.000. Sheet 3, Ariège. Geological Institute, Leiden University. Leidsche Geologische Mededelingen, 22, 351-418.
De Sitter, L.U., Zwart, H.J., Kleinsmiede, W.F.J., 1960. Geological map of the Central Pyrenees, 1/50.000. Sheet 4, Valle de Aran. Geological Institute, Leiden University.
Delgado, J., 1993. Caracterización mineralógica, físico-química y geoquímica de los skarns del contacto norte del batolito de la Maladeta (Val d’Aran, Lleida). Doctoral Thesis. Universitat de Barcelona, 412pp.
Druschel, G.K., Rosenberg, P.E., 2001. Non-magmatic fracture- controlled hydrothermal systems in the Idaho Batholith: South Fork Payette geothermal system. Chemical Geology, 173, 271-291.
Dunne, T., 1978. Water in environmental planning. New York (USA), W.H. Freeman and Company, 818pp.
ENHER, 1985. Estudio geológico y geoquímico de detalle de la zona de surgencias termales de Caldes de Bohí y Artiés. Permiso de exploración nº 3952 “Pallars”. Barcelona, 1985.
Fernández, M., Banda, E., 1989. An approach to the thermal field in northeastern Spain. Tectonophysics, 164, 259-266.
Fournier, R.O., 1979. Geochemical and hydrological considerations and the use of enthalpy-chloride diagrams in the prediction of underground conditions in hotspring systems. Journal of Volcanology and Geothermal Research, 5, 1-16.
Gimeno, M.J., Auqué, L.F., Gómez, J.B., García, I., 2007. Geochemical features and geothermometric modeling of Les Escaldes geothermal system (Andorra). In: Bullen, T., Wang, Y. (eds.). Proceedings of the 12th International Symposium on water-rock interaction, 321-324.
Gimeno, M.J., Auqué, L.F., Acero, P., Gómez, J.B., 2014. Hydrogeochemical characterisation and modelling of groundwaters in a potential geological repository for spent nuclear fuel in crystalline rocks (Laxemar, Sweden). Applied Geochemistry, 45, 50-71.
Gómez, J.B., Gimeno, M.J., Auqué, L.F., Acero, P., 2014. Characterisation and modelling of mixing processes in groundwaters of a potential geological repository for nuclear wastes in crystalline rocks of Sweden. Science of the Total Environment, 468-469, 791-803.
Gupta, S.C., Wang, D., Dong, H.K., 2006. Water Retention in Soil. In: Lal, R. (ed.). Encyclopedia of Soil Science. Boca Ratón, Florida, CRC Press, 2nd Edition, 1864-1869.
Gutenbrunner, C., Bender, T., Cantista, P., Karagülle, Z., 2010. A proposal for a worldwide definition of health resort medicine, balneology, medical hydrology and climatology. International Journal of Biometeorology, 54, 495-507.
Han, D.M., Liang, X., Jin, M.G., Currell, M.J., Song, X.F., Liu, C.M., 2010. Evaluation of groundwater hydrochemical characteristics and mixing behavior in the Daying and Qicun geothermal systems, Xinzhou Basin. Journal of Volcanology and Geothermal Research, 189, 92-104.
IGME, 1984. Proyecto de investigación geotérmica preliminar del Pirineo Central, zona meridional del Prelitoral Catalán e Islas Baleares. Unpublished, Vol. 3, Ministerio de Industria y Energía, 188pp.
Jones, A., Montanarella, L., Jones, R. (eds.), 2005. Soil Atlas of Europe. Luxembourg, Office for Official Publications of the European Communities, 128pp.
Kele, S., Breitenbach, S.F.M., Capezzuoli, E., Meckler, A.N., Ziegler, M., Millan, I.M., Kluge, T., Deák, J., Hanselmann, K., John, C.M., Yan, H., Liu, Z., Bernasconi, S.M., 2015. Temperature dependence of oxygen- and clumped isotope fractionation in carbonates: a study of travertines and tufas in the 6–95ºC temperature range. Geochimica et Cosmochimica Acta, 168, 172-192.
Kleinsmiede, W.F.J., 1960. Geology of the Valle de Arán (Central Pyrenees). Leidsche Geologische Mededelingen, 25, 129-245.
Lal, R., Shukla, M. (eds.), 2004. Principles of Soil Physics. New York and Basel, Marcel Dekker Inc., 682pp.
López-Moreno, J.L., García-Ruiz, J.M., 2004. Influence of snow accumulation and snowmelt on stream flow in the central Spanish Pyrenees. Hydrological Sciences Journal, 49(5), 787-802.
Martinez, J.L., Raiber, M., Cendónc, D.I., 2017. Using 3D geological modelling and geochemical mixing models to characterise alluvial aquifer recharge sources in the upper Condamine River catchment, Queensland, Australia. Science of the Total Environment, 574, 1-18.
Mey, P.H.W., 1965. Geological map of the Ribagorzana and Baliera valleys, Central Pyrenees, 1/25.000. Geological Institute, Leiden University.
Mey, P.H.W., 1968. Geology of the Upper Ribagorzana and Tor Valleys, sheet 8, Central Pyrenees, Spain. Leidsche Geologische Mededelingen, 41, 229-292.
Michard, G., 1983. Recueil de données thermodynamiques concernant les équilibres eaux-minéraux dans les réservoirs géothermaux. Rapport EUR 8590 FR, 156 p.
Michard, G., Fouillac, C., 1980. Contrôle de la composition chimique des eaux thermals sulfurées sodiques du Sud de la France. In: Tardy, Y. (ed.). Geochimique des interactions entre les eaux, les minéraux et les roches. Elements Tarbes, 147-166.
Michard, G., Roekens, E., 1983. Modelling of the chemical composition of alkaline hot waters. Geothermics, 12, 161-169.
Michard, G., Sanjuan, B., Criaud, A., Fouillac, C., Pentcheva, E.N., Petrov, P.S., Alexieva, R., 1986. Equilibria and geothermometry in hot waters from granites of SW Bulgaria. Geochemical Journal, 20, 159-171.
Moreno, L., Azcón, A., Navarrete, P., Garrido, E., 1997. Caracterización geológica e hidroquímica de las surgencias con influencia termal en los macizos graníticos de la zona axial pirenaica. Soria, VII Congreso de Geoquímica de España, 392-399.
Ninot, J.M., Ferré, A., 2008. Plant diversity across five vegetation belts in the Pyrenees (Catalonia, Spain). Collectanea Botanica (Barcelona), 27, 65-74.
Oyarzún, R., Barrera, F., Salazar, P., Maturana, H., Oyarzún, J., Aguirre, E., Alvarez, P., Jourde, H., Kretschmer, N., 2014. Multimethod assessment of connectivity between surface water and shallow groundwater: the case of Limarí River basin, north-central Chile. Hydrogeology Journal, 22, 1857-1873.
Parkhurst, D.L., Appelo, C.A.J., 1999. User’s Guide to PHREEQC (Version 2), a computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. U.S. Geological Survey, Water- Resources Investigations Report, 99-4259.
Sanjuan, B., Michard, A., Michard, G., 1988. Influence of the temperature of CO2-rich springs on their Al and REE contents. Chemical Geology, 68, 57-67.
Savage, D., 2011. A review of analogues of alkaline alteration with regard to long-term barrier performance. Mineralogical Magazine, 75, 2401-2418.
Soulé, J.C., 1990. Circulations profondes en milieu granitique: eaux sulfurées des Pyrénées. Hydrogeologie, 4, 297-299.
Thornthwaite, C.W., 1948. An approach toward a rational classification of climate. Geographical Review, 38, 55-94
Todorovic, B.Z., Stojilkovic, D.T., Petrovic Pantic, T., Mitic, N.C., Nikolic, L.S., Cakic, S., 2016. Hydrochemistry and aragonite scaling in the Sijarinska spa (Serbia). Carbonates and Evaporites, 31, 367-374.
Van Manen, S.M., Wallin, E., 2012. Ground temperature profiles and thermal rock properties at Wairakei, New Zealand. Renew Energy, 43, 313-321.
Van Middlesworth, P.E., Wood, S.A., 1998. The aqueous geochemistry of the rare earth elements and yttrium. Part 7. REE, Th and U contents in thermal springs associated with the Idaho Batholith. Applied Geochemistry, 13, 861-884.
Vilaplana, J.M., 1983. Quaternary glacial geology of Alta Ribargorça basin (Central Southern Pyrenees). Acta Geológica Hispánica, 18(3-4), 217-233.
Vinchon, C., 1977. Contribution à l’étude pétrographique du Silurien des Pyrénées centrales espagnoles (région du Rio Esera, Province de Huesca et région de Llavorsi, Province de Lerida). D.E.A., Lille University, 172pp.
Zandvliet, J., 1960. The geology of the Upper Salat and Pallaresa valleys, Central Pyrenees, France/Spain 1/20.000. Leidsche Geologische Mededelingen, 25, 1-127.
Zwart, H.J., 1977. Six cross sections through the Central Pyrenees, 1/50.000. Geological Institute, Leiden University.
Zwart, H.J., 1979. The Geology of the Central Pyrenees. Leidsche Geologische Mededelingen, 50, 1-74.
Zwart, H.J., Roberti, K.F., 1976. Geological map of the Central Pyrenees. 1:50.000. Sheet 9, Flamisell-Pallaresa. Geological Institute, Leiden University.
Michard, G., 1990. Behaviour of major elements and some trace elements (Li, Rb, Cs, Sr, Fe, Mn, W, F) in deep hot waters from granitic areas. Chemical Geology, 89, 117-134.
ELECTRONIC APPENDIX I



n.a.: not analyzed ; Total alkalinity is expressed as HCO3

n.a.: not analyzed; b.l.d.: below detection limit; Total alkalinity is expressed as HCO3 -

(*)2003 data; n.a.: not analyzed; b.l.d.: below detection limit; Total alkalinity is expressed as HCO3

n.a.: not analyzed; b.l.d.: below detection limit; Total alkalinity is expressed as HCO3

(*)2003 data; n.a.: not analyzed; b.l.d.: below detection limit; Total alkalinity is expressed as HCO3 -