Major-and Trace-Element Geochemistry of Geothermal Water from the Nappe Zone, Northern Tunisia: Implications for Mineral Prospecting and Health Risk Assessment

: A comprehensive hydrogeochemical survey of the geothermal waters from the Nappe Zone (Maghrebides fold-and-thrust belt) was undertaken to determine the origins of geothermal waters and to assess the health risks associated with their potentially toxic elements. A total of 11 geothermal water and 3 stream water samples were collected and analysed for major and trace elements (As, B, Ba, Fe, Mn, Pb, Sr, Zn). Two main geothermal water groups were highlighted by hydrogeochemical diagrams and multivariate analyses (PCA, HCA): the ﬁrst group is the Na − Cl type, TDS > 10 g/L, controlled by deep circulation, while the second group is the Na-Cl − HCO 3 type, TDS < 2 g/L, and controlled by shallow circulation. A curved hydrogeochemical evolution path, observed from mixed bicarbonate shallow groundwater to chloride geothermal water, indicates that the interaction with evaporites drives the chemistry of the geothermal samples. On these, the As enrichments come from sulphide oxidation polymetallic mineralisation during the upwelling to the surface from E–W major lineaments. Therefore, E–W lineaments are potential areas for mineral prospecting. The health risk assessment reveals that the concentration of potentially toxic elements in geothermal waters are lower than the guideline values for the protection of freshwater aquatic life and dermal exposure (bathing or balneology).


Introduction
Geothermal and mineral waters have been a worldwide interest in recent years [1,2].The importance given to these natural resources is mainly related to their suitability for various purposes such as bathing, therapy, heating and cooling, and irrigation.Tunisia is known to have several potential geothermal sites, which are a result of its location on the African Plate and the Mediterranean Ridge.
On a global scale, drinking waters distributed for human consumption through the water supply network exhibit a Ca-HCO 3 chemistry and Total Dissolved Solids (TDS) < 1 g/L (freshwater), characteristic of shallower groundwater that dissolves limestone, which is the most abundant rock constituting aquifers [3].However, within the study area, drinking waters show a dominant Na-Cl chemistry and an average TDS value of 1234 ± 154 mg/L (brackish salinity) [4].This could be due to the presence of saline deposits and/or mixing with deeper, higher-salinity waters as geothermal waters.
In Tunisia, geothermal waters are mainly employed for balneotherapy and the irrigation of crops in greenhouses.Several researchers have inventoried different geothermal water springs [5][6][7][8], evaluated their potential [9,10], described their chemical and isotopic compositions [1,11], and investigated their relation with the geological and structural patterns of Tunisia [12,13].Previous studies also revealed the predominance of crustal gas signatures, the dissolution of halite-bearing Triassic evaporites, and the influence of regional geology, in particular the structural control of the faults, which allow for the deep circulation and rising of fluids on geothermal springs across Tunisia [1,10,[12][13][14][15].
The occurrence of potentially toxic elements (PTEs) in geothermal waters represents a serious environmental concern due to the high toxicity, long persistence, and bioaccumulation potential of such elements [16,17].The most common PTEs in geothermal waters come from geogenic inputs, including soil leaching, weathering processes, hydrothermal activity and ore deposits, and anthropogenic inputs, including domestic and industrial effluents and mining wastes [18][19][20].These PTEs can enter the human body mainly through direct ingestion and dermal contact during bathing and balneology [21,22].The discharge of geothermal waters in the environment may also pose a serious threat to human health and ecosystems through the contamination of aquatic ecosystems (groundwater and stream water) [20,23].Therefore, in other countries, some studies based on the hydrochemistry of geothermal water have been conducted for environmental surveys rather than for mineral exploration purposes [24,25].
Though the traditional geochemistry methods taking stream sediments and soils as mediums are less useful for deep, concealed ore-deposit prospecting, the hydrochemistry of spring water can rapidly identify deep-buried mineralisation in a large region given an excellent vertical migration ability and limited budget.Geothermal waters ascend from a geothermal reservoir and emerge at the surface, providing reliable information about the subsurface environments.Several previously published manuscripts have described the use of hydrochemistry for mineral prospecting in ore deposits [19,[26][27][28].
Since knowledge of the occurrence of trace elements in geothermal water is still limited in Tunisia, a comprehensive hydrogeochemical survey of the geothermal waters from the Nappe Zone, northern Tunisia, was undertaken by the Tunisian National Office of Mines (ONM).In this study, the physical parameters (pH, temperature T, and electrical conductivity EC) and major (Na + , K + , Ca 2+ , Mg 2+ , HCO 3 − , Cl − , NO 3 − , and SO 4 2− ) and trace elements (As, B, Ba, Fe, Mn, Pb, Sb, Sr, and Zn) in 11 geothermal water and 3 stream water samples were measured.The Nappe zone has high geothermal potential due to its location on the Maghrebides fold-and-thrust belt.This region is already known to host a variety of Pb-Zn (Ba-Sr-F-Fe-Hg) hydrothermal deposits and is assumed to have some of the highest subsurface mineral potential in Tunisia.This study aims to apply hydrochemistry and multivariate statistics through hierarchical cluster analysis (HCA) and principal component analysis (PCA) of geothermal water for mineral prospecting and health risk assessment.

Geological, Hydrogeological, Geothermal, and Mineral Settings
The Nappe Zone of northern Tunisia is part of the Maghrebide-Alpine belt (Figure 1), stretching from the Moroccan Rif in the west to southern Italy to the east [29,30].The Maghrebide-Alpine belt has its genesis linked to the collision between the African and Eurasian plates, which occurred during the Langhian [31,32].The Nappe Zone exhibits a heterogeneous geology represented by Miocene folds, and the Tellian thrust sheets over the Atlasic foreland to the south.The Tellian thrust sheets are mainly composed by upper Cretaceous-Eocene shales and limestones [33,34], and are overlain by the Numidian thrust sheet, which is made up of Oligocene-Burdigalian claystones and sandstones [34].The thrust domain tectonically overhangs the Atlasic foreland, consisting of limestone, dolomite, and clay.The Atlasic foreland is represented by large-scale NE-SW synclines and anticlines, which are crossed by SE graben and horst.The NE-SW and E-W trending faults have facilitated the ascension and effusion of saliferous Triassic formations consisting mainly of gypsum, halite, and anhydrite [34][35][36].From the Middle to Late Miocene, magmatism favoured the emplacements of felsic and mafic bodies in the Nefza, Tabarka, and Mogod regions [33,37].Consequent to the Tortonian NW-trending compressive event, Neogene post-nappes basins start to develop and were filled by lacustrine limestones and volcanoclastic deposits [38].
Northern Tunisia is characterised by deep basins and faults formed by a series of extensional and compressional tectonics that have reactivated inherited structures from the Cretaceous [35,36].The main deep faults and lineaments are the NE-SW and NW-SE trending faults and a series of N 80 • E lineaments (Figure 1) [39,40].Moreover, some shallow faults and lineaments are recognised consequent to Neogene thrusting [41].
The intense faulting has triggered the deep circulation of hydrothermal fluids and the formation of a wide variety of ore deposits [38,[42][43][44][45]. Based on a geochemical survey of stream sediments in the Nappe Zone, it was determined that the Mississippi Valley-type Pb-Zn ore deposits are primarily associated with NE-SW trending faults, whereas the Pb-As-Sb sedimentary exhalative (SEDEX) type deposits are predominantly controlled by E-W-trending faults within the Neogene continental basins [46].As and Sb are attributed to the hydrothermal activity triggered by the Late Miocene magmatism [39].Arsenopyrite, realgar, and orpiment are common ore minerals associated with SEDEX ore deposits [39,47].High As concentrations are also found in galena (0.1-0.6%) from the Jalta ore deposits (Nappe Zone) [48].Therefore, As is an useful pathfinder for the exploration of Pb-Zn SEDEX ore deposits in the Nappe Zone [46].
In hydrogeological terms, the Miocene-Quaternary continental deposits covering the accumulation basins favour the infiltration of meteoric waters which percolate to the deeper levels (from the Lias to the Quaternary), where tectonic activity plays an important role [49].The intense tectonic activity provides ideal conditions for the formation of numerous geothermal systems in the Nappe Zone: deep circulation and upward movement of fluids through basins characterised by an average geothermal gradient of 35 • C km −1 [7,8,50] mainly resulted from the superposition of multiple geological units reaching more than 3000 m thick [34].According to previous studies [1,7,8], spring waters are mainly meteoric in origin, with the dominance of the sodium-chloride type.
Finally, according to the Köppen-Geiger climate classification system, the studied area is classified as a Warm Mediterranean Climate (Csa) with a mean annual air temperature of 18.

Sample Collection and Analytical Procedure
A total of 11 geothermal water samples, homogenously distributed over the Nappe Zone, were collected in March 2022.In addition, 3 stream water samples were also collected from the main tributaries of Oued Medjerda, largely draining the study region.The water samples were collected in double pre-cleaned bottles (polypropylene), one acidified (0.5 mL HNO3 65% Suprapur ® Supelco in 100 mL sample) for major cation and trace element analyses and the other stored without acidification for anion analyses.Physical parameters including temperature (T), pH, and electrical conductivity (EC) were measured in the field using a multiparameter probe (Hanna Instruments, Smithfield, RI, USA).Geothermal fluid samples were filtered through 0.45 µm Millipore filters and the concentrations of major cations, calcium Ca 2+ , magnesium Mg 2+ , sodium Na + , and potassium K + were determined with a Flame Atomic Absorption Spectrometer (PinAAcle 500, PerkinElmer, Waltham, MA, USA).The anion concentrations of chloride (Cl − ) and bicarbonate (HCO3 − ) were obtained using titration methods.Nitrate (NO3 − ) and sulphate (SO4 2− ) were obtained via gravimeter method.Total dissolved solids (TDS) parameter was obtained via calculation method [54,55].Standard and blank solutions were prepared and repetitively run after ten samples.The concentrations of As, B, Ba, Fe, Mn, Pb, Sr, and Zn were measured using inductively coupled plasma atomic emission spectroscopy (Ultima C HORIBA, Jobin Yvon GmbH, Oberursel, Germany).The detection limits for trace elements ranged from 1 to 5 µg L −1 and from 0.1 to 1 mg L −1 for major elements.[29,39,43,53], and location of the sampling sites.

Sample Collection and Analytical Procedure
A total of 11 geothermal water samples, homogenously distributed over the Nappe Zone, were collected in March 2022.In addition, 3 stream water samples were also collected from the main tributaries of Oued Medjerda, largely draining the study region.The water samples were collected in double pre-cleaned bottles (polypropylene), one acidified (0.5 mL HNO 3 65% Suprapur ® Supelco in 100 mL sample) for major cation and trace element analyses and the other stored without acidification for anion analyses.Physical parameters including temperature (T), pH, and electrical conductivity (EC) were measured in the field using a multiparameter probe (Hanna Instruments, Smithfield, RI, USA).Geothermal fluid samples were filtered through 0.45 µm Millipore filters and the concentrations of major cations, calcium Ca 2+ , magnesium Mg 2+ , sodium Na + , and potassium K + were determined with a Flame Atomic Absorption Spectrometer (PinAAcle 500, PerkinElmer, Waltham, MA, USA).The anion concentrations of chloride (Cl − ) and bicarbonate (HCO 3 − ) were obtained using titration methods.Nitrate (NO 3 − ) and sulphate (SO 4 2− ) were obtained via gravimeter method.Total dissolved solids (TDS) parameter was obtained via calculation method [54,55].Standard and blank solutions were prepared and repetitively run after ten samples.The concentrations of As, B, Ba, Fe, Mn, Pb, Sr, and Zn were measured using inductively coupled plasma atomic emission spectroscopy (Ultima C HORIBA, Jobin Yvon GmbH, Oberursel, Germany).The detection limits for trace elements ranged from 1 to 5 µg L −1 and from 0.1 to 1 mg L −1 for major elements.
The quality of the analytical procedure was verified based on analysing replicate samples (10%) and reference water samples (from the National Office of Mines, Tunis, Tunisia).
According to the standard method 1030E, checking analyses' correctness [54], the accuracy of the major element concentrations was checked by calculating the ionic-balanceerror (IBE) between total concentration of cations and anions, as follows: Finally, the recovery result of trace elements (Concentration measured/Concentration certified × 100) ranged from 89 to 107% depending on the element considered.

Multivariate Statistical Analysis and Geochemical Codes
One important objective of the present study was to identify the origin apportionment of the measured hydrochemical variables in the study region.Multivariate statistical analysis tools, such as hierarchical cluster analysis (HCA) and principal component analysis (PCA), are effective techniques for the exploration of hydrochemical datasets [56][57][58].In the current study, PCA and HCA were applied to 17 variables (pH, EC, TDS, Na + , K + , Ca 2+ , Mg 2+ , HCO 3 − , Cl − , NO 3 − , SO 4 2− , B, Ba, Mn, Pb, Sr, and Zn).Arsenic and Fe were excluded from the multivariate statistical analysis because their concentrations were below the detection limit in more than 50% of the samples.The centred-log-ratio (clr) transformation was performed to the raw data prior to multivariate statistical analysis in order to overcome the problem of the compositional nature of hydrochemical data [59].HCA and PCA were computed with R 3.5.1 (R Foundation for Statistical Computing, Vienna, Austria) using rgr [60], factoextra [61], and FactoMineR [62] packages.
Hierarchical cluster analysis (HCA) is a common multivariate technique widely used to classify variables or observations within a hydrochemical dataset based upon their similarities [56,58].A matrix of distance measures or similarity between objects (variables or observations) is firstly computed and then the most similar objects are linked to produce new clustered objects.In this study, HCA was carried out with data of a matrix based on Ward's method of linkage and Euclidean distance between observations as the measure of similarity [63], which assumed to produce the most distinctive hydrochemical clusters [64].
Principal component analysis (PCA) was successfully used to reduce numerous hydrochemical variables into few principal components (PCs) based on a correlation or covariance matrix [56,57].The geochemical interpretation of these principal components provides insight into the major processes that can control the distribution of hydrochemical variables.
Regarding the choice of geochemical software, mineral saturation indexes and CO 2 fugacity were calculated for each sample using two different codes: The Geochemist's Workbench ® suite of programs, release 12 (hereafter GWB code) [65], and the USGS's PHREEQC-Interactive Windows Interface, version 3 (hereafter, PhreeqcI code) [66].This approach allows for an estimation of the error, as there may be differences due to the varying methods of calculating species activity in a dissolution reaction (Q), as well as differences in equilibrium constants.Therefore, for activity calculations, both the B-dot and the Pitzer model were considered [67], in combination with different thermodynamic databases: thermo.tdatand thermo_ymp.R2.tdat with GWB code, and llnl.dat and pitzer.datwith PhreeqcI code.The GWB code was also used to classify water samples on the basis of equivalent concentration of major dissolved constituents, whereas PHREEQCI code was used to model the sulphate-water interaction.Finally, the temperature at depth of geothermal water samples was inferred by the multicomponent geothermometry computer program GeoT, version 2.1 [68].

Geothermal Water Quality Assessment
Health risk assessment models proposed by USEPA are widely used to evaluate the non-carcinogenic and carcinogenic risks associated with potentially toxic elements in aquatic ecosystems [69].Potential non-carcinogenic risks through ingestion and dermal contact are assessed based on Hazard Index (HI) and hazard quotient (HQ) using the following equations: Hazard Index(HI) = HQ ing + HQ der Hazard Quotient (HQ) ing = CDI ing /RfD ing Hazard Quotient (HQ) der = CDI der /RfD der where CDI is chronic daily intake (mg kg −1 day −1 ) and C 0 is the concentration of trace elements in a water sample (µg kg −1 ).The exposure parameters are given in Table 1, while the values of oral reference dose (RfD; mg kg −1 day −1 ), dermal permeability constant (K p ; cm h −1 ), and Gastrointestinal absorption (GIABS; unitless) are shown in Table 2.
In this study, the RfD and K p values were taken from [69,70], respectively.As no RfD value was proposed for Pb, the provisional levels provided by [71] were used.
HQ and HI values lower than 1 indicate no significant non-cancer health risk, while values of HQ and HI greater than 1 depict existing likelihood of non-cancer health effects occurring, and the probability increases as the values rise [69,72].

Results
The results of the physicochemical analyses are reported in Table 3. Subsequently, the findings will be discussed using various approaches.

Hydrochemical Features
The hierarchical cluster analysis, Langelier-Ludwig's [76] and Spencer's [77,78] diagrams are suitable tools to understand and identify hydrogeochemical features, types, and evolution of water.Moreover, the following descriptors were also used in this study: (i) "water type" or "geochemical facies" denomination furnished by the GWB code and using the highest cation and anion concentration in an equivalent basis of a sample; (ii) TDS salinity classification [79]; and (iii) temperature classification [80].The HCA of the geothermal and stream water samples identified three main clusters at the phenon line with a linkage distance of about two (Figure 2): cluster C1, including three samples, which can be divided into two sub-clusters C1.1 (which includes only the sample #1 Sidi Abdelkader) and C1.2 (#3 H. Sayala and #10 H. S. Ali Daoua); cluster C2, containing four samples, and cluster C3, which can be divided into two sub-clusters C3.1 and C3.2, containing three samples (stream water) and four samples, respectively.

Hydrochemical Features
The hierarchical cluster analysis, Langelier-Ludwig s [76] and Spencer s [77,78] diagrams are suitable tools to understand and identify hydrogeochemical features, types, and evolution of water.Moreover, the following descriptors were also used in this study: (i) "water type" or "geochemical facies" denomination furnished by the GWB code and using the highest cation and anion concentration in an equivalent basis of a sample; (ii) TDS salinity classification [79]; and (iii) temperature classification [80].The HCA of the geothermal and stream water samples identified three main clusters at the phenon line with a linkage distance of about two (Figure 2): cluster C1, including three samples, which can be divided into two sub-clusters C1.1 (which includes only the sample #1 Sidi Abdelkader) and C1.2 (#3 H. Sayala and #10 H. S. Ali Daoua); cluster C2, containing four samples, and cluster C3, which can be divided into two sub-clusters C3.1 and C3.2, containing three samples (stream water) and four samples, respectively.The thermal water samples from the first cluster (C1) are characterised by the highest values of Na + , Ca 2+ , Mg 2+ , K + , Cl − , SO4 2− , and TDS (Tables 3 and 4) with Cl-Na (Na/Cl < 1) type saline to brine waters (Figure 3a).The thermal water samples from the first cluster (C1) are characterised by the highest values of Na + , Ca 2+ , Mg 2+ , K + , Cl − , SO 4 2− , and TDS (Tables 3 and 4) with Cl-Na (Na/Cl < 1) type saline to brine waters (Figure 3a).In particular, the brine sample #10, H. S. Ali Daoua, is the saltiest (TDS = 58.4g/L) and has a Na-Cl composition (Na/Cl > 1).The brackish samples from cluster C2 show (i) a weakest value of Ca 2+ , Mg 2+ , Na + , Cl − , and TDS and (ii) a main hydrochemical type is Na−Cl with a trend towards the Na-HCO 3 field, represented by sample #11 (Figure 3a).The brackish samples of cluster C3 have the highest concentrations of HCO 3 − , and the main hydrochemical Na-Cl types persist if the dominant cations and anions are looked at singularly.However, in the sub-cluster C3.1, represented by hypothermal stream water samples, samples could be also alternatively classified as Ca(Mg)-Cl, because the sum Ca + Mg in eq/L is higher than the Na + K concentration.In a similar way, sub-cluster C3.2 shows a Na-Cl/Ca-Cl composition with a lower Ca(Mg)-Cl component, except for fresh sample #2 which is a fresh, Ca-HCO 3 , orthothermal water suggesting a trend towards the typical shallow groundwater field (Figure 3a).Geothermal water can mix with shallow groundwater or infiltrating water during their upward movement through deep faults and fractures, producing bicarbonate-rich waters [81].Cluster C3, containing stream water samples, and cluster C2 are related to the weakest linkage distance (around three), suggesting a higher affinity of mixing, while both these clusters (C3 and C2) join with C1 at a higher linkage distance (6), indicating a lower hydrochemical affinity.
The geochemical signature of different clusters can be further discriminated using the Langelier-Ludwig classification diagram, where the equivalent percentages of (Ca + Mg), (Na + K), (HCO 3 ), and (Cl + SO 4 ) are projected.Consistent with the findings of HCA, the initial clusters are clearly defined in terms of the distribution of the samples on the Langelier-Ludwig diagram (Figure 3a) and according to the distribution of the proposed genetic water fields on that diagram [82].  1, as described in the text, and previously published water samples on the same sites.In (A): elliptic field depicts possible water origins [82]; dotted path depicts the sequential dissolution of evaporate (first gypsum and then halite) [77].In (B): hatched field corresponds to brine waters, which dissolve an evaporitic assemblage formed by 1.0 CaSO4 and 0.5 CaSO4 [77].In both: coloured and open symbols are from this study and the literature [1,8], respectively; S-hexagon depicts the composition of Mediterranean seawater (SW) [83].[77,78] (B) diagrams (equivalent basis) of geothermal and stream water samples classified according to hierarchical cluster analysis.The numbers on the symbols refer to the specific samples in Table 1, as described in the text, and previously published water samples on the same sites.In (A): elliptic field depicts possible water origins [82]; dotted path depicts the sequential dissolution of evaporate (first gypsum and then halite) [77].In (B): hatched field corresponds to brine waters, which dissolve an evaporitic assemblage formed by 1.0 CaSO 4 and 0.5 CaSO 4 [77].In both: coloured and open symbols are from this study and the literature [1,8], respectively; S-hexagon depicts the composition of Mediterranean seawater (SW) [83].
Samples from C2 and C3, characterised by various water facies and brackish salinity, are largely distributed in the mixing zone of Figure 1a, reflecting the interaction with shallow groundwater and saline/brine water.C1.2 samples fall in the upper left corner of the Langelier-Ludwig diagram, suggesting that these highly mineralised Na−Cl thermal waters are mainly controlled by the halite dissolution.Sample C1.1 is represented by Na + −Cl − rich thermal water from the coastal region of Lake Ichkeul (Figure 1).The hydro- dynamics of the lake are influenced by input/output contributions from both freshwater streams and saltwater from Bizerte lagoon [84].Despite this, previous investigations of the isotope's composition revealed that C1.1 had a meteoric signature [1].Therefore, the position of the C1.1 sample within the "marine water" field in Figure 3a could be random.Instead, the composition could be alternatively explained as the dissolution of gypsum and, subsequently, halite by a meteoric water [77].Indeed, both samples of this study and from previous published studies [1,8] fall on the curved evolutive path describing that process (Figure 3a).In fact, the composition of sample C1.1, as well as the other C1.2 samples in Figure 3b, is similar to brines resulting from evaporite dissolution, as also demonstrated by the undersaturation in gypsum and anhydrite (Table S1).For the Sayala sample, there could be contributions from other sulphate minerals (e.g., glauberite), as observed in other Mediterranean brines [77].
In terms of trace element contents, Ba, Mn, Pb, Sr, and Zn concentrations of the geothermal waters in the study region are higher than those of stream waters (Table 3).This could be related to the higher reactivity and ability of geothermal water to dissolve trace elements from the host rocks during circulation [2,85,86].In a different way, the trace element concentrations of stream waters may be controlled by the continuous dilution of surface waters with meteoric waters.The geothermal samples in sub-cluster C1.2 are the only two samples having an As concentration of 140 and 339 µg/L, which is higher than the guideline value of 10 µg/L [87].

Hydrogeochemical Processes
Geochemical ratios are widely carried out to characterise the main processes and mechanisms that control water geochemistry [88,89].The plot of Na + versus Cl − reveals that geothermal water samples in all clusters are categorised by a molar Na + /Cl − ratio close to one, indicating that the dissolution of evaporitic formations rich in halite dominates the dissolved Na + and Cl − in the geothermal waters (Figure 4a).This finding is consistent with previous studies, revealing that the chemistry of Tunisian geothermal water is primarily governed by the dissolution of the saliferous Triassic formations rich in halite [1,8,15].
The correlation of SO 4 2− vs. Ca 2+ shows that the majority of samples are below and along the dissolution line of gypsum and/or anhydride (Figure 4b), reflecting the important contribution of both of the minerals to the geothermal water chemistry in the Nappe Zone [1].Samples from sub-cluster C1.2 are distributed far from this line, with SO 4 2− /Ca 2+ ratios higher than two.This deviation from the gypsum dissolution line may be explained by either (i) a decrease in Ca 2+ contents due to the ion exchange process involving silicate weathering; (ii) an increase in sulphate contents due to sulphide oxidation [2,58]; or (iii) the dissolution of sulphate minerals different from gypsum and anhydrite [77].The HCO 3 − vs. (Ca 2+ + Mg 2+ ) correlation (Figure 4c) reveals sub-cluster samples C1.2 are plotted below the dissolution line of the carbonates (calcite and dolomite), with HCO 3 − /(Ca 2+ + Mg 2+ ) ratios lower than one, indicating an increase in Ca 2+ and Mg 2+ contents.Therefore, the dissolution of Ca(Mg)-bearing sulphate minerals is plausible in those C1.2 samples.In a different way, samples from cluster C2 are distributed beneath the dissolution line of the carbonates, confirming an ion exchange process commonly resulting from silicate weathering.These findings are supported by the (Ca 2+ + Mg 2+ )-(HCO 3 − + SO 4 2− ) vs. (Na + + K + )-Cl − plot.Indeed, in that diagram (Figure 4d), C2 samples are plotted along the line, characteristic of Ca 2+ + Mg 2+ deficiency and Na enrichment, confirming the ion exchange process.Coastal geothermal water sample (C1.1) plots along the opposite side of this line as a response to (i) a slight seawater contribution or (ii) the reverse ion exchange process triggered by the marine intrusion, as already observed in several coastal aquifers in Tunisia [88,90].Samples from sub-cluster C1.2 are distributed far from that line.In particular, the large Na excess in brine sample #10 could be due to the dissolution of glauberite, Na 2 Ca(SO 4 ) 2 .
Environments 2023, 10, 151 12 of 26 or (ii) the reverse ion exchange process triggered by the marine intrusion, as already observed in several coastal aquifers in Tunisia [88,90].Samples from sub-cluster C1.2 are distributed far from that line.In particular, the large Na excess in brine sample #10 could be due to the dissolution of glauberite, Na2Ca(SO4)2.Chloride and boron are conservative elements that are slowly affected by secondary processes, including participation in weathering minerals, and then are regarded as indicators providing evidence of mixing [91,92].In this study, the Cl − /B ratio of the geothermal waters in clusters fluctuated between 2000 and 10,000 (Figure 4e).These high ratios support the fact that water-rock interactions are the main process in the evolution of water chemical composition [91,92].In addition, a high positive correlation between B and Cl concentrations (r = 0.81) of the geothermal waters from C1.1, C3, and C3 was observed (Figure 4e), indicating mixing between the geothermal and stream waters.
A plot of the NO 3 − /Cl − ratio vs. Cl − is commonly used to identify the main origins of NO 3 − [88].Figure 4f shows that the samples from C2 and C3 exhibit the highest NO 3 − /Cl − ratio and the lowest Cl concentrations, suggestion that NO 3 − could mainly result from nitrate fertilisers.Groundwater contamination by the agricultural inputs is widely reported in Tunisia [88].This result is another support for the mixing process.
The clustering results includenitially identified from HCA (Figure 2) are projected on the PC1-PC2 plane (Figure 5), upon which their chemical signature will be defined.The initial clusters are once again obviously defined in terms of the distribution of the water samples (Figure 5).In terms of PC1, its main effect is the opposite between saline and fresh waters.Cluster 1 displays positive loadings on PC1, and it corresponds to the most saline samples, poor in HCO 3 with a low pH, which supports the idea that C1 is mainly controlled by the dissolution of halite for C1.2 and the seawater intrusion for the sample collected near the coast (C1.1).These results are also in full agreement with a previously published study [1].Conversely, the samples from clusters 2 and 3 reveal negative loadings on PC1, and they contain abundant HCO 3 , Pb, Zn, Mn, and Ba with a high pH.This may be the result of mixing and the dilution of geothermal fluids with shallow groundwater rich in HCO 3 − , as previously revealed by HCA and ionic ratio results.The high contents of HCO 3 − may reflect the effect of carbonate formations, largely occurring in northern Tunisia, while the increase in trace elements (Pb, Zn, Mn, and Ba) may derive from the weathering of Pb-Zn ore deposits and mining wastes [93][94][95].These trace elements are the main ore-forming elements associated with Pb-Zn deposits hosted by carbonates (limestone and dolomites) and Miocene shales in the Nappe Zone [39,46].
The antithetical correlation between K + and Ca 2+ − Mg 2+ through PC2 is attributed to the different genesis of the waters.Sub-cluster C1.2 and cluster 2 occur on the positive axis of PC2 with abundant K + , indicating that evaporite dissolution and ion exchange processes play important roles in the control of the solutes in geothermal water for sub-cluster C1.2 and cluster 2, respectively.
Sub-cluster C3.1 (stream water samples) shows positive loadings on PC3 (Figure S1), and this corresponds to samples that are abundant in B, while C3.2 occurs on the positive side corresponding to samples with abundant Sr and a high electrical conductivity (EC).Therefore, boron enrichment may result from the flushing of clayey materials with freshwater, causing bivalent cation depletion [96,97].The cation exchange process is accompanied by a decrease in ionic strength [97], which is consistent with the results of the current study.
PC4 distinguishes sub-cluster C3.1 (stream water samples), rich in Sr, from C3.2 (geothermal samples), with abundant NO3 − (Figure S2).This may occur as C3.1 denotes natural stream water rich in Sr with a lower NO3 − , while C3.2 represents shallower geothermal water contaminated by NO3 − , supporting the evidence found from the plot NO3 − /Cl − vs. Cl − (Figure 4f).This component can be attributed to the effect of mixing and dilution processes.In fact, increasing the mixing fraction induces the decreasing of nitrate concentration via the dilution effect.

Geothermometry
The geothermometric evaluation of waters dissolving evaporites (in this case, salt deposits of seawater origin mainly of the Triassic age) is challenging, because most of the chemically based geothermometers published so far assume equilibrium with silicate minerals.In the geothermal squared diagram of Figure 6a [98], it can be observed that the waters from cluster C1, which exhibit surface temperatures in line with the classic "hot water" meaning (>36 °C) [99], are closer to the full equilibrium curve than C2.Specifically, according to [8], it is noticeable that the temperature at depth does not exceed 100 °C.This is further confirmed by the application of the Na-K-Ca geothermometer [100,101], which is based on silicate equilibrium, and where the R factor > 50 indicates a temperature not significantly different from the surface temperature.Subtracting the measured PC4 distinguishes sub-cluster C3.1 (stream water samples), rich in Sr, from C3.2 (geothermal samples), with abundant NO 3 − (Figure S2).This may occur as C3.1 denotes natural stream water rich in Sr with a lower NO 3 − , while C3.2 represents shallower geothermal water contaminated by NO 3 − , supporting the evidence found from the plot NO 3 − /Cl − vs. Cl − (Figure 4f).This component can be attributed to the effect of mixing and dilution processes.In fact, increasing the mixing fraction induces the decreasing of nitrate concentration via the dilution effect.

Geothermometry
The geothermometric evaluation of waters dissolving evaporites (in this case, salt deposits of seawater origin mainly of the Triassic age) is challenging, because most of the chemically based geothermometers published so far assume equilibrium with silicate minerals.In the geothermal squared diagram of Figure 6a [98], it can be observed that the waters from cluster C1, which exhibit surface temperatures in line with the classic "hot water" meaning (>36 • C) [99], are closer to the full equilibrium curve than C2.Specifically, according to [8], it is noticeable that the temperature at depth does not exceed 100 • C.This is further confirmed by the application of the Na-K-Ca geothermometer [100,101], which is based on silicate equilibrium, and where the R factor > 50 indicates a temperature not significantly different from the surface temperature.Subtracting the measured temperature values to those obtained by the Na-K-Ca geothermometer, a range of 60-100 • C is obtained (Table S3).The only outlier is represented by brine sample #10, which shows the highest deep temperature of 116 • C.However, this latter value is not plausible due to the significant contribution of cations from evaporites to this water.Indeed, using the simultaneous equilibrium of anhydrite, gypsum, and calcite, a temperature of 64 • C is obtained using the GeoT code [68] (Table S3).Such a discrepancy is due to the dissolution of evaporites, whose dissolved cations can lead to distorted geothermometric information when applied to a silicate geothermometer as Na-K-Ca.Using the PHREEQCI software, version 3 [66], to model the water-evaporite dissolution, it is possible to check the plausibility of the temperatures as depicted in Figure 6a.The SO 4 2− /Cl − versus SO 4 plot is very useful to trace sulphate origin and its fate in groundwater [102,103].In Figure 6b, brackish to saline/brine samples increase sulphate concentration with decreasing SO 4 2− /Cl − ratios, which is the opposite behaviour of the traced gypsum/anhydrite dissolution curves.However, it should also be noted that sulphate dissolution could be enhanced by high ionic strength furnished by sodium and chloride (i.e., from halite dissolution), but also more soluble sulphate minerals such as glauberite.Indeed, saline to brine samples of the C1 cluster fall on or close to the different interaction path traced using the PHREEQCI code, the deep inferred temperature between 40 and 100 • C, and increasing Na + = Cl − mol concentration.
Environments 2023, 10, 151 15 of 26 temperature values to those obtained by the Na-K-Ca geothermometer, a range of 60-100 °C is obtained (Table S3).The only outlier is represented by brine sample #10, which shows the highest deep temperature of 116 °C.However, this latter value is not plausible due to the significant contribution of cations from evaporites to this water.Indeed, using the simultaneous equilibrium of anhydrite, gypsum, and calcite, a temperature of 64 °C is obtained using the GeoT code [68] (Table S3).Such a discrepancy is due to the dissolution of evaporites, whose dissolved cations can lead to distorted geothermometric information when applied to a silicate geothermometer as Na-K-Ca.Using the PHREEQCI software, version 3 [66], to model the water-evaporite dissolution, it is possible to check the plausibility of the temperatures as depicted in Figure 6a.The SO4 2− /Cl − versus SO4 plot is very useful to trace sulphate origin and its fate in groundwater [102,103].In Figure 6b, brackish to saline/brine samples increase sulphate concentration with decreasing SO4 2− /Cl − ratios, which is the opposite behaviour of the traced gypsum/anhydrite dissolution curves.However, it should also be noted that sulphate dissolution could be enhanced by high ionic strength furnished by sodium and chloride (i.e., from halite dissolution), but also more soluble sulphate minerals such as glauberite.Indeed, saline to brine samples of the C1 cluster fall on or close to the different interaction path traced using the PHREEQCI code, the deep inferred temperature between 40 and 100 °C, and increasing Na + = Cl − mol concentration.  .In (A), AC-pentagon depicts the isochemical dissolution of average crust [98], SW, the Mediterranean seawater, as in Figure 3.The mixing between AC and SW or full equilibrium water at different temperatures are also shown (dashed).In (B), bicarbonate waters are from [98,102,103].Sulphate dissolution modelling by PHREEQCI code [66] is as shown curves with different designs.In particular, from left to right: gypsum dissolution at 18 °C and Na + = Cl − = 2.8 mmol/l (dash double-dot), gypsum dissolution at 40 °C and Na + = Cl − =26 mmol/L (dash), anhydrite dissolution at 100 °C and Na + = Cl − = 200 mmol (dash dot), anhydrite0.98-glauberite0.02dissolution at 80 °C and Na + = Cl − = 880 mmol (dot).In both, symbols as in Figure 3. Open triangles are samples from previously published studies [1,8].

Conceptual Model of Geothermal Fluid and Its Implication for Mineral Prospection
By integrating both the results of this study and those of previous surveys [1,8], a conceptual model of the geothermal system in the Nappe Zone is illustrated in the schematic diagram (Figure 7).The geothermal waters are primarily of meteoric origin.When meteoric waters infiltrate from the high latitudes through faults, waters are conductively heated and react with evaporite-carbonate rocks, forming dominant Na-Cltype water.However, due to the difference in tectonic patterns, there are two major groups of geothermal systems in the study region; one is controlled by deep circulation (Group I), and the other is controlled by shallow circulation (Group II).
The first group (Group I) includes geothermal water samples from cluster C1.2 that are characterised by the highest values of TDS, Na + , Ca 2+ , Mg 2+ , Cl − , and SO4 2− .Additionally, these springs are located close to the E-W major lineaments, providing channels for the migrations of infiltrating meteoric water to the deep crust.High TDS values of these geothermal waters can be explained by the long residence times, extensive water-rock interactions, long flow pathways, and limited mixing with shallow

Conceptual Model of Geothermal Fluid and Its Implication for Mineral Prospection
By integrating both the results of this study and those of previous surveys [1,8], a conceptual model of the geothermal system in the Nappe Zone is illustrated in the schematic diagram (Figure 7).The geothermal waters are primarily of meteoric origin.When meteoric waters infiltrate from the high latitudes through faults, waters are conductively heated and react with evaporite-carbonate rocks, forming dominant Na-Cl-type water.However, due to the difference in tectonic patterns, there are two major groups of geothermal systems in the study region; one is controlled by deep circulation (Group I), and the other is controlled by shallow circulation (Group II).flushing (dilution), which induced the precipitation of As [110].The amount of As in wate can be also retained in the soils and stream sediments via adsorption reactions into iro oxyhydroxides/oxides [110], similar to the widespread As anomalies observed in stream sediments in north Tunisia [46].This could explain the low Fe concentrations i geothermal water [113].

Geothermal Water Quality Assessment
Compared with the water quality guidelines for the protection of freshwater aquati life, trace-element concentrations in most geothermal springs are within the criterio maximum concentrations (CMC) set by USEPA [22] (Table 5).The concentrations o different trace elements in all geothermal samples are less than the guideline values fo irrigation water [114], with the exception of As and Sr from two and all springs respectively.When compared with water quality guidelines for drinking water, trace element concentrations in most geothermal springs exceed the drinking water guidelin values established by the WHO [87,115] and the European Union [116,117].
Currently, the geothermal water rich in potentially toxic elements in the study are is mainly used by local populations for external purposes such as spas and balneology which may induce significant adverse effects on human health by way of skin contact.Th assessment results of non-carcinogenic health risks to children and adults exposed to trac elements in geothermal waters via dermal contact are shown in Figure 8a,b.HQ exceeding one (HQs ≥ 1) indicate harmful health effects, while HQs lower than one (HQ < 1) suggest a tolerable hazard level [21].Overall, the HQDermal values estimated fo adults are above those for children (Figure 8a,b).Furthermore, the HQDermal for differen geothermal water samples indicate that both groups of the population are exposed to PTE in the following order: As > Mn > Pb > Sr > Ba > Zn, where As in geothermal water sample from sub-cluster C1.2 (H.Sidi Ali Daoua and H. Sayala) is the element with the highes HQDermal values in both populations.HQDermal values are found to be below than on The first group (Group I) includes geothermal water samples from cluster C1.2 that are characterised by the highest values of TDS, Na + , Ca 2+ , Mg 2+ , Cl − , and SO 4 2− .Additionally, these springs are located close to the E-W major lineaments, providing channels for the migrations of infiltrating meteoric water to the deep crust.High TDS values of these geothermal waters can be explained by the long residence times, extensive water-rock interactions, long flow pathways, and limited mixing with shallow groundwater.Therefore, these spring waters could be evolved in separate and isolated geothermal systems.Based on a geophysical approach, a previous study revealed that geothermal systems at H. Sayala spring are confined by a Miocene clay layer [104].When the meteoric water migrates downward along E-W major lineaments, it should dissolve anhydrite and gypsum imbedded along the faults and/or in deeper host rocks, similar to results found in north-eastern Algeria [85].The estimated reservoir temperature in Sidi Ali Daoua geothermal spring is between 100 and 120 • C. Assuming an average geothermal gradient of 35 • C km −1 [7,8,50], it is possible that a deep aquifer is located at a depth of ~3 km below the surface, where geothermal waters reach a thermo-chemical equilibrium with the halite and gypsum rich evaporites.Such depth and temperature estimations are consistent with those estimated for geothermal waters from northern Tunisia [15] and Algeria [85].Deep-seated faults reaching a 2000-3000 m depth were already identified in the Nappe Zone by Essid et al.'s studies [35,36].This group (Group I) is also typified by As enrichments, which are probably related to long residence times and high interactions of the deep geothermal water with the wall-rocks [105,106].The H. Sidi Ali Daoua and H. Sayala As-rich geothermal waters coincide with Oued Maaden and Sidi Bou Aouane hydrothermal ore deposits, respectively, where As is found as a major component in sulphide minerals such as arsenopyrite (FeAsS), realgar (AsS), and orpiment (As 2 S 3 ) and as a secondary component in sulphur minerals such as chalcopyrite (CuFeS 2 ) and pyrite (FeS) [37,40,45].In northern Tunisia, many Pb-As-Sb sedimentary exhalative (SEDEX) ore deposits are distributed along E-W major lineaments and are hosted by Miocene shales in continental basins [46].These basins are furthermore characterised by EW-trending planar As-Sb anomalies resulting from weathering and the mobilisation of As-bearing ore bodies [39,46].Therefore, the regional E-W major lineaments have served as ore-conducting pathways and ore-hosting networks.
The upwelling of deep geothermal water rich in dissolved arsenic and sulphate along these faults could be related to the oxidation of sulphide minerals hosted in Miocene shales constituting the cap-rocks of these deep geothermal springs.During the rising of geothermal water from the deeper reservoir to the surface, As is affected by a series of biotic and abiotic processes such as oxidation, reduction, and As-S redox cycling.Arsenic oxidation can occur under both aerobic and anaerobic environments via bacterial processes [19,92].High sulphate concentrations in the As-rich geothermal waters support the oxidation of As-bearing minerals and the resulting acidity is buffered by the dissolution of carbonates, which is supported by the high Ca 2+ concentrations, particularly for H. Sidi Ali Daoua (1150 mg/L).In fact, according to numerous previous studies [107,108] the neutral-alkaline pHs of deep spring waters are explained by the fact that deeper-circulating waters are more likely to react with carbonate layers and, thus, be buffered.Additionally, it was shown that the neutral spring waters are more prone to leach arsenic from the host rocks than acid spring waters [109,110].In contrast, cationic PTEs (Pb, Zn and Fe), commonly associated with ore deposits, display relatively low concentrations in the As-rich geothermal waters.This can be explained by the fact that the high carbonate contents of sedimentary deposits in northern Tunisia buffer the solution at a neutral-alkaline pH, causing the immobilisation of cationic PTEs [111,112].Therefore, in the case of a neutral-alkaline environment like northern Tunisia, arsenic, sulphate, and calcium anomalies detected in geothermal waters are likely to be indicators for mineralised zones.
The second group (Group II) includes exclusively relatively diluted mixed springs (TDS < 1500 mg/ L) across the Nappe Zone (clusters C2 and C3), as discussed earlier.Thus, these spring waters are mainly linked to shallower faults and fractures rather than major and deep lineaments (Figure 7).In fact, according to [8], the reservoir depth is estimated at 1 km below the ground surface, while the reservoir temperature is about 65 • C for H. Bourguiba spring.The low TDS values and lack of high As concentrations in these deeply mixed spring waters explain the short-term water-rock exchanges at relatively moderate temperatures.In a similar way, a lack of elevated As concentrations was described in the mixed spring water in Yellowstone National Park, Wyoming, USA, by the effect of massive flushing (dilution), which induced the precipitation of As [110].The amount of As in water can be also retained in the soils and stream sediments via adsorption reactions into iron oxyhydroxides/oxides [110], similar to the widespread As anomalies observed in stream sediments in north Tunisia [46].This could explain the low Fe concentrations in geothermal water [113].

Geothermal Water Quality Assessment
Compared with the water quality guidelines for the protection of freshwater aquatic life, trace-element concentrations in most geothermal springs are within the criterion maximum concentrations (CMC) set by USEPA [22] (Table 5).The concentrations of different trace elements in all geothermal samples are less than the guideline values for irrigation water [114], with the exception of As and Sr from two and all springs, respectively.When compared with water quality guidelines for drinking water, trace-element concentrations in most geothermal springs exceed the drinking water guideline values established by the WHO [87,115] and the European Union [116,117].
Currently, the geothermal water rich in potentially toxic elements in the study area is mainly used by local populations for external purposes such as spas and balneology, which may induce significant adverse effects on human health by way of skin contact.The assessment results of non-carcinogenic health risks to children and adults exposed to trace elements in geothermal waters via dermal contact are shown in Figure 8a,b.HQs exceeding one (HQs ≥ 1) indicate harmful health effects, while HQs lower than one (HQs < 1) suggest a tolerable hazard level [21].Overall, the HQDermal values estimated for adults are above those for children (Figure 8a,b).Furthermore, the HQDermal for different geothermal water samples indicate that both groups of the population are exposed to PTEs in the following order: As > Mn > Pb > Sr > Ba > Zn, where As in geothermal water samples from sub-cluster C1.2 (H.Sidi Ali Daoua and H. Sayala) is the element with the highest HQDermal values in both populations.HQDermal values are found to be below than one in all springs (Figure 8a,b), which indicates that geothermal water in this study region is safe for human health in terms of external use. in all springs (Figure 8a,b), which indicates that geothermal water in this study region is safe for human health in terms of external use.However, changing dietary customs, especially under changing climate scenarios such as water scarcity, increasing temperature, and diversity in rainfall patterns, can incite the local populations to consume these geothermal waters, placing their health at risk.Arsenic in geothermal waters from H. Sidi Ali Daoua and H. Sayala springs is the element with the highest HQIngestion values in both populations (Figure 8c,d), while Pb HQIngestion for children in all geothermal water samples is above one.As a result, geothermal waters are unsuitable for drinking and both of the populations could experience health risks and significant hazard levels (Figure 8c,d).The HI values (sum of HQ values) of the selected trace elements are mainly greater than one for both groups of the population and are mostly driven by the ingestion of Pb and As, particularly in H. Sidi Ali Daoua and H. Sayala springs (Figure 9).HI values for children are greater than those for adults, reflecting that children are more sensitive than adults when exposed to geothermal water rich in PTEs.However, changing dietary customs, especially under changing climate scenarios such as water scarcity, increasing temperature, and diversity in rainfall patterns, can incite the local populations to consume these geothermal waters, placing their health at risk.Arsenic in geothermal waters from H. Sidi Ali Daoua and H. Sayala springs is the element with the highest HQIngestion values in both populations (Figure 8c,d), while Pb HQIngestion for children in all geothermal water samples is above one.As a result, geothermal waters are unsuitable for drinking and both of the populations could experience health risks and significant hazard levels (Figure 8c,d).The HI values (sum of HQ values) of the selected trace elements are mainly greater than one for both groups of the population and are mostly driven by the ingestion of Pb and As, particularly in H. Sidi Ali Daoua and H. Sayala springs (Figure 9).HI values for children are greater than those for adults, reflecting that children are more sensitive than adults when exposed to geothermal water rich in PTEs.[116]; f Food and Agricultural Organization (FAO) for irrigation water [114]; -no limit recommended.

Conclusions
Geothermal water in the Nappe Zone is fed by meteoric water infiltrating into the subsurface, where it is conductively heated from the surrounding rocks at an anomalous geothermal gradient (35 °C/Km) due to the complex structural and geological settings.In this context, infiltrating meteoric water has been evolved in two main different geothermal systems: deep and superficial.
Water springs from superficial systems having chemical compositions similar to those of local stream water with relatively low concentrations of TDS (<2 g/L) appear to be related to short circulation pathways and mixing processes with shallow groundwater

Conclusions
Geothermal water in the Nappe Zone is fed by meteoric water infiltrating into the subsurface, where it is conductively heated from the surrounding rocks at an anomalous geothermal gradient (35 • C/Km) due to the complex structural and geological settings.In this context, infiltrating meteoric water has been evolved in two main different geothermal systems: deep and superficial.
Water springs from superficial systems having chemical compositions similar to those of local stream water with relatively low concentrations of TDS (<2 g/L) appear to be related to short circulation pathways and mixing processes with shallow groundwater often contaminated by trace elements and nitrate.
Water springs from deep and confined systems having high concentrations of TDS (>10 g/L) appear to be related to deeply circulated meteoric waters that have been severely modified by long water-rock interactions at temperatures ranging not higher than 100 • C. The compositions of these spring waters suggest they are fed by separate and isolated geothermal systems which are controlled by deep faults.The regional E-W major lineaments have provided pathways for the infiltration of meteoric water into the deep crust (more than 3000 m).The deepest geothermal waters at H. Sidi Ali Daoua and H. Sayala springs show the highest arsenic, sulphate, and calcium concentrations, which are probably evidence of bacterially mediated arsenic oxidation and reduction.Arsenic is commonly associated with Pb-Zn SEDEX ore deposits mainly hosted by Miocene shales that constitute the cap-rocks of the deepest geothermal springs.Therefore, the Pb-As-Sb SEDEX deposit type was successfully detected by the hydrochemistry of deep confined geothermal water.
This study evidences that in the case of a neutral-alkaline environment like northern Tunisia, As, SO 4 2− , and Ca 2+ can be regarded as useful pathfinders for exploring the Pb-Zn SEDEX ore deposits in the Nappe Zone.Arsenic shows a higher mobility in neutralalkaline pH environments, consistent with competitive sorption with HCO 3 − .In contrast, cationic PTEs (e.g., Pb, Zn, Fe) are immobilised in the alkaline pH range due to the high carbonate content.

Supplementary Materials:
The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/environments10090151/s1,Table S1: Mineral saturation indexes and CO 2 fugacity results obtained using GWB and PHREEQCI codes and different thermodynamic databases; Table S2: Component loadings of four principal components, eigenvalues, and variances; Figure S1: PCA loadings of the hydrochemical parameters on the PC1-PC3, and PCA scores of water samples classified according to hierarchical cluster analysis; Figure S2: PCA loadings of the hydrochemical parameters on the PC1-PC4 and PCA scores of water samples classified according to hierarchical cluster analysis; Table S3: Application of Na-K-Ca Fournier and Truesdell's geothermometer [100,101], and comparison with GeoT multimineral code [68].

Figure 2 .
Figure 2. Main clusters in the HCA s dendrogram for the geothermal and stream water samples.The different colours highlight main clusters and sub-clusters.The black line is "phenon line".

Figure 2 .
Figure 2. Main clusters in the HCA's dendrogram for the geothermal and stream water samples.The different colours highlight main clusters and sub-clusters.The black line is "phenon line".

Figure 3 .
Figure 3. Langelier-Ludwig s[76] (A) and Spencer s[77,78] (B) diagrams (equivalent basis) of geothermal and stream water samples classified according to hierarchical cluster analysis.The numbers on the symbols refer to the specific samples in Table1, as described in the text, and previously published water samples on the same sites.In (A): elliptic field depicts possible water origins[82]; dotted path depicts the sequential dissolution of evaporate (first gypsum and then halite)[77].In (B): hatched field corresponds to brine waters, which dissolve an evaporitic assemblage formed by 1.0 CaSO4 and 0.5 CaSO4[77].In both: coloured and open symbols are from this study and the literature[1,8], respectively; S-hexagon depicts the composition of Mediterranean seawater (SW)[83].

Figure 3 .
Figure 3. Langelier-Ludwig's[76] (A) and Spencer's[77,78] (B) diagrams (equivalent basis) of geothermal and stream water samples classified according to hierarchical cluster analysis.The numbers on the symbols refer to the specific samples in Table1, as described in the text, and previously published water samples on the same sites.In (A): elliptic field depicts possible water origins[82]; dotted path depicts the sequential dissolution of evaporate (first gypsum and then halite)[77].In (B): hatched field corresponds to brine waters, which dissolve an evaporitic assemblage formed by 1.0 CaSO 4 and 0.5 CaSO 4[77].In both: coloured and open symbols are from this study and the literature[1,8], respectively; S-hexagon depicts the composition of Mediterranean seawater (SW)[83].

Figure 5 .
Figure 5. PCA loadings of the hydrochemical parameters on the first two principal components, and PCA scores of water samples classified according to hierarchical cluster analysis.

Figure 5 .
Figure 5. PCA loadings of the hydrochemical parameters on the first two principal components, and PCA scores of water samples classified according to hierarchical cluster analysis.

Figure 8 .
Figure 8. Hazard quotient (HQ) for dermal contact (a,b) and ingestion (c,d) in children and adults, respectively.

Figure 8 .
Figure 8. Hazard quotient (HQ) for dermal contact (a,b) and ingestion (c,d) in children and adults, respectively.

Figure 9 .
Figure 9. Hazard Index (HI) for children and adults.

Table 1 .
Exposure parameters used for the health risk assessment.

Table 3 .
Chemical analysis of geothermal and stream water samples classified according to hierarchical cluster analysis (major elements and TDS concentrations in mg/L; trace elements concentrations, from As to Zn, are in µg/L; Nd: not detected).

Table 5 .
Permissible limits for drinking water quality and aquatic life protection (µg/L).

Table 5 .
Permissible limits for drinking water quality and aquatic life protection (µg/L).