The Evolution and Sources of Major Ions in Hot Springs in the Triassic Carbonates of Chongqing, China

: Thermal groundwater in the Lower and Middle Triassic carbonates in Chongqing, China, is mainly concentrated in anticlines. Hot springs (32.9 to 57 ◦ C) with SO 4 -Ca type waters and Total Dissolved Solids (TDS) of 1620 to 2929 mg / L emerge in the middle and the plunging ends of the structures. Multivariate methods are used to analyze the hydrochemical characteristics of the waters, and identify the sources of the main dissolved components, providing an insight into the evolution of the environment in which they formed. Hierarchical cluster analysis of compositional data di ﬀ erentiates samples in the study area into three categories: high TDS-high Ca 2 + and SO 42 − water; medium TDS-high Na + and Cl − water; and low TDS-high HCO 3 − water. Factor analysis and ion ratio relationships show that Ca 2 + and SO 42 − are mainly derived from the dissolution of gypsum and anhydrite within the geothermal reservoir, with some addition of SO 42 − from coal-bearing cap rocks. The main source of HCO 3 − , is in the dissolution of dolomite and CO 2 that also promotes the incongruent dissolution of albite and K-feldspar, adding Na + and K + to the groundwater. Reverse modelling of the transfers of each phase shows, in three models, that the minerals dissolved decrease progressively—with the exception of halite and albite. Combined with the hydrochemical characteristics of hot water in the same reservoir in the adjacent area (Cl-Na type, TDS of 13.37 g / L), a process of desalination of the hot water can be conﬁrmed, which has not yet reached the ‘freshwater’ stage dominated by HCO 3 − .


Introduction
Subsurface brines, saline springs, and hot waters (32-64 • C) are common in the Triassic carbonates of the Sichuan Basin in southwestern China [1][2][3], in the eastern and northeastern parts of the basin [3,4]. The hot springs of the eastern basin, in the Chongqing area, appear in the axes, cores, two flanks and plunging ends of a series of northeast-southwest (NE-SW) trending anticlines ( Figure 1). Geothermal wells have been drilled near the hot springs and on the flanks of the anticlines, and hot water is also present in the tunnels in these areas. The thermal groundwaters are characteristically of SO 4 -Ca type, with Total Dissolved Solids (TDS) of 2 to 3 g/L [4,5]. The utilization and development of the thermal waters in this area has grown rapidly in the past two decades, and thus an understanding of the hydrochemical characteristics and compositional evolution of the geothermal water-apparently strongly controlled by the anticlines-is important in providing insight into the geothermal potential of other carbonate reservoirs.  Much research has been applied to the thermal waters of Chongqing, using correlation analysis and traditional graphical methods (Piper, Schoeller, Ludwig-Langelier and Durov diagrams). Xiao et al. (2013) discussed the water-rock interactions of the Bei hot spring in the reservoir using δ 13 C-HCO 3 − and δ 34 S-SO 4 2− . These reveal that the main water rock reaction is the dissolution of gypsum, followed by the reaction of CO 2 with the host rock to produce HCO 3 − [6]. The hot water in Chongqing reservoirs flows along the flanks of the anticlines, obtaining heat from below and recharging Water 2020, 12, 1194 3 of 20 in the cores of the anticlines from infiltration of precipitation (rain) in mountain areas, at approximately 670 to 1500 m altitude, where the carbonates outcrop [4]. The thermal groundwater issues at the surface as hot springs in the axis or plunging ends of the anticlines [5,6]. However, the circulation compositional evolution and the quantitative analysis of water-rock interactions have been little studied. Methods widely applied in the analysis of the hydrogeochemical evolution of geothermal waters include major ion ratios, isotope tracers, graphic methods (Gibbs diagrams, box diagrams, etc.), multivariate statistical analysis and hydrogeochemical simulations. Many researchers have studied the evolution characteristics and water-rock interactions of hot waters at local or regional scales using isotope tracers of carbon (δ 13 C-CO 2 , δ 13 C-TDIC) [7,8], beryllium (δ 11 B) [9,10], sulphur and strontium (δ 34 S-SO 4 2− , 87 Sr/ 86 Sr) [11,12]. Gibbs (1970) proposed a descriptive model (a Gibbs diagram) based on the relationship between TDS and Na + /(Na + + Ca 2+ ), and Cl − /(Cl − + HCO 3 − ) to explain the mechanism controlling the source of dissolved components in groundwater, from atmospheric precipitation, rock weathering, processes of evaporation, and crystallization process [13]. This model has become the classical model for the study of hydrochemistry and has been widely used worldwide [14][15][16].
Because of their advantages in processing large-scale hydrochemical data sets, multivariate statistical methods have been successfully applied in the study of hydrogeochemical processes [17][18][19][20]. They have become powerful tools for characterizing the spatial and temporal variations of hydrogeochemistry and identifying the sources of chemical components. In this paper, R-model factor analysis (FA) is used to study the relationship between variables, and Q-model cluster analysis (CA), to group samples. Based on the principle of mass conservation, hydrogeochemical simulation can quantitatively calculate the amount of material transfer during water-rock interaction and the dissolution and precipitation state of various minerals, mapping the evolution of the groundwater [21][22][23]. The graphics of multimineral stability are used to analyze the main water-rock interaction in a groundwater system from the perspective of mineralogy and its application in hydrogeochemical research can provide more intuitive evidence [24][25][26]. The purpose of this study is to discuss the unique hydrochemical characteristics of the thermal groundwaters of Chongqing and the compositional evolution of the major ions in the geothermal reservoir, using field and historical data and combining the main hydrochemical methods with simulation models. The results provide an understanding of water-rock interactions in these Triassic carbonate reservoirs and may serve as a reference for similar sites elsewhere.

Description of the Study Area
The Sichuan Basin, lies in the northwestern part of the Upper Yangtze Platform, with extensional faulting defining the approximately rhombic basin margins ( Figure 1) [27]. During the Palaeozoic, the basin received wide-spread marine deposits, including limestones. During the Mesozoic, these peripheral faults were transformed by compression, and the basin was inverted, and deposits were predominantly lacustrine. Cenozoic deposits consisted principally of intermountain plain fluvial facies [28].
The study area in the eastern part of the basin lies in a high steep fold belt, formed by the early Cretaceous-Eocene Sichuan Movement. A series of NE-NNE trending and nearly parallel anticlinal folds were developed in this area with broad hills forming parallel ridge-valleys ( Figure 1). Mesozoic rocks form the main bodies of the folds, characterized by long narrow axes, high cores, and steep flanks. Middle and Lower Triassic carbonates outcrop in the axes. By contrast, complementary synclines are wide and gentle and in these areas the Triassic rocks are deeply buried. Some of the anticlines are asymmetrical with dips on adjacent limbs varying from relatively gentle 20-30 • , to steeper 40-70 • , with some even vertical or inverted [1]. The Wentangxia anticline ( Figure 1) extends more than 200 km NE-NNE along its axis. The northern part of the anticline is steep and narrow with a width of 5 to 8 km and an elevation of 800 to 110 m. Limestones with typical karst valleys are exposed along the axis. In the Guanyinxia anticline, dips are gentler to the east than to the west. This anticline is rich in geothermal water and wells have been drilled along both flanks. The Tongluoxia anticline is Water 2020, 12, 1194 4 of 20 slightly curved, with an axis of approximately 230 km, with the southern end plunging obliquely and connecting with the eastern flank of the Nanwenquan anticline [29]. The Yangtze River, forming the base level of erosion and discharge in the study area, cuts across the anticlines to form deep gorges and, together with the tributary Jialingjiang River, and forms a deep gorge with drainage points for thermal groundwater (up-flow springs and erosion springs) [1,6].
The geothermal reservoirs in this area are mainly in carbonate rocks of the Lower Triassic Jialingjiang Group (T 1 j) and Middle Triassic Leikoupo Group (T 2 l), together about 600 m thick, and outcrop in the cores of the anticlines (Figure 2). The Jianglingjiang Group can be divided into four sections (T 1 j 1 , T 1 j 2 , T 1 j 3 , and T 1 j 4 ). T 1 j 1 consists of grey limestone, and mudstone interbedded with limestone bearing insect tracks and bio-limestone. T 1 j 2 , comprises grey argillaceous dolomite with limestone (the lower part of the T 1 j 2 is gypsum), grey-green and purple mudstone. T 1 j 3 , consists of light grey medium to thick bedded limestone, argillaceous limestone, dolomitic limestone and bio-limestone. T 1 j 4 , comprises light grey dolomite, argillaceous dolomite, and limestone. The overlying Leikoupo Group (T 2 l) is mainly limestone and yellowish-grey dolomite, with argillaceous dolomite, limestone with anhydrite (the bottom part of the T 2 l is hydromica clay rock). The Upper Triassic Xujiahe Group (T 3 xj) comprises siliciclastic rocks (450 m thick) and outcrops on the flanks of the anticlines. These include grey-white feldspar-rich quartz arenite, quartz arenite with grey-black shale, carbonaceous shale and a thin coal seam. The impervious shale and coal at the bottom of T 3 xj and red Jurassic sandy mudstones (1000 m or more in thickness) outcrop on both sides of the anticlines and constitute the cap rocks of the geothermal reservoirs ( Figure 3). The Lower Triassic Feixianguan Group (T 1 f) that lies beneath the reservoirs consists of purple-red, grey-purple and yellow-green mudstones and dark-grey marl with varying thickness. This acts as an impermeable base, and locally outcrops in the cores of some of the anticlines (e.g.; the Wentangxia and Guanyinxia anticlines, Figure 1). deep gorges and, together with the tributary Jialingjiang River, and forms a deep gorge with drainage points for thermal groundwater (up-flow springs and erosion springs) [1,6]. The geothermal reservoirs in this area are mainly in carbonate rocks of the Lower Triassic Jialingjiang Group (T1j) and Middle Triassic Leikoupo Group (T2l), together about 600 m thick, and outcrop in the cores of the anticlines (Figure 2). The Jianglingjiang Group can be divided into four sections (T1j 1 , T1j 2 , T1j 3 , and T1j 4 ). T1j 1 consists of grey limestone, and mudstone interbedded with limestone bearing insect tracks and bio-limestone. T1j 2 , comprises grey argillaceous dolomite with limestone (the lower part of the T1j 2 is gypsum), grey-green and purple mudstone. T1j 3 , consists of light grey medium to thick bedded limestone, argillaceous limestone, dolomitic limestone and biolimestone. T1j 4 , comprises light grey dolomite, argillaceous dolomite, and limestone. The overlying Leikoupo Group (T2l) is mainly limestone and yellowish-grey dolomite, with argillaceous dolomite, limestone with anhydrite (the bottom part of the T2l is hydromica clay rock). The Upper Triassic Xujiahe Group (T3xj) comprises siliciclastic rocks (450 m thick) and outcrops on the flanks of the anticlines. These include grey-white feldspar-rich quartz arenite, quartz arenite with grey-black shale, carbonaceous shale and a thin coal seam. The impervious shale and coal at the bottom of T3xj and red Jurassic sandy mudstones (1000 m or more in thickness) outcrop on both sides of the anticlines and constitute the cap rocks of the geothermal reservoirs ( Figure 3). The Lower Triassic Feixianguan Group (T1f) that lies beneath the reservoirs consists of purple-red, grey-purple and yellow-green mudstones and dark-grey marl with varying thickness. This acts as an impermeable base, and locally outcrops in the cores of some of the anticlines (e.g.; the Wentangxia and Guanyinxia anticlines, Figure 1).

Sampling and Testing
A total of 28 hot water samples were analysed in this study. Of these, 18 were sampled in the field survey in October 2013 and August 2016, while analyses of the other 10, cited in [30], were collected from geothermal wells in the main urban area of Chongqing in the wet season of 2009 for comparison. Variable hydrochemical parameters, including temperature (±0.1 °C) and pH (±0.1 unit), were measured in-situ with a hand-held thermometer (TP101), pH and ORP meters (CT-6821), calibrated in the field prior to sampling. Water samples were stored in polyethylene bottles, including one large (1500 mL) bottle and two smaller (550 mL) bottles that were rinsed with water samples three times before actual sampling. SiO2 was measured at the laboratory of the Beijing Brigade of Hydrogeology and Engineering Geology by Si-Mo yellow colorimetry (UV spectrophotometer). Laboratory testing of other chemical parameters was completed at the Analysis Center of the Beijing Research Institute of Uranium Geology, determining the following parameters: concentrations of cations K + , Na + , Ca 2+ and Mg 2+ , and anions Cland SO4 2− were determined using ion chromatography (ICS-1100 and 883 Basic IC plus); HCO3 − by the volumetric method (HCl titration using an AT-510 automatic titration analyzer); and Li + , and Sr 2+ using inductively coupled plasma-mass spectrometry (NexION300D). The results of anion-cation equilibrium tests show that the equilibrium errors of the water samples ranged from 2.92% to 3.63%.

Cluster Analysis (CA)
Clustering analysis is an unsupervised pattern recognition technology that can reveal the intrinsic structure or potential features of data sets without any prior assumption. It classifies the objects of the system into categories or clusters based on their proximity or similarity, and is commonly used in the management of hydrogeochemical data [31][32][33]. Q-mode hierarchical cluster analysis (HCA) is the most widely used method. In this, comparisons based on multiple parameters from different samples are made and the samples are grouped in clusters according to their "similarity" to each other, forming higher order clusters step-by-step [34][35][36]. The Euclidean distance may be used quantitatively to describe the similarity between samples. In other words, the difference between the statistical values of samples can be expressed by "distance". Without considering other factors, groundwaters with a close "distance" should have the same recharge source or closely hydraulic connection. The further the "distance" is, the greater the difference in chemical composition and the weaker the hydraulic connection becomes [37,38]. HCA was performed using Ward's linkage rule. This method uses an analysis of variance approach to evaluate the distance between clusters, attempting to minimize the sum of squares of any two (hypothetical) clusters that could be formed at each step [39,40]. The results of HCA are finally shown by a dendrogram, which intuitively reflects

Sampling and Testing
A total of 28 hot water samples were analysed in this study. Of these, 18 were sampled in the field survey in October 2013 and August 2016, while analyses of the other 10, cited in [30], were collected from geothermal wells in the main urban area of Chongqing in the wet season of 2009 for comparison. Variable hydrochemical parameters, including temperature (±0.1 • C) and pH (±0.1 unit), were measured in-situ with a hand-held thermometer (TP101), pH and ORP meters (CT-6821), calibrated in the field prior to sampling. Water samples were stored in polyethylene bottles, including one large (1500 mL) bottle and two smaller (550 mL) bottles that were rinsed with water samples three times before actual sampling. SiO 2 was measured at the laboratory of the Beijing Brigade of Hydrogeology and Engineering Geology by Si-Mo yellow colorimetry (UV spectrophotometer). Laboratory testing of other chemical parameters was completed at the Analysis Center of the Beijing Research Institute of Uranium Geology, determining the following parameters: concentrations of cations K + , Na + , Ca 2+ and Mg 2+ , and anions Cland SO 4 2− were determined using ion chromatography (ICS-1100 and 883 Basic IC plus); HCO 3 − by the volumetric method (HCl titration using an AT-510 automatic titration analyzer); and Li + , and Sr 2+ using inductively coupled plasma-mass spectrometry (NexION300D). The results of anion-cation equilibrium tests show that the equilibrium errors of the water samples ranged from 2.92% to 3.63%.

Cluster Analysis (CA)
Clustering analysis is an unsupervised pattern recognition technology that can reveal the intrinsic structure or potential features of data sets without any prior assumption. It classifies the objects of the system into categories or clusters based on their proximity or similarity, and is commonly used in the management of hydrogeochemical data [31][32][33]. Q-mode hierarchical cluster analysis (HCA) is the most widely used method. In this, comparisons based on multiple parameters from different samples are made and the samples are grouped in clusters according to their "similarity" to each other, forming higher order clusters step-by-step [34][35][36]. The Euclidean distance may be used quantitatively to describe the similarity between samples. In other words, the difference between the statistical values of samples can be expressed by "distance". Without considering other factors, groundwaters with a close "distance" should have the same recharge source or closely hydraulic connection. The further the "distance" is, the greater the difference in chemical composition and the weaker the hydraulic connection becomes [37,38]. HCA was performed using Ward's linkage rule. This method uses an analysis of variance approach to evaluate the distance between clusters, attempting to minimize the sum of squares of any two (hypothetical) clusters that could be formed at each step [39,40]. The results of HCA are finally shown by a dendrogram, which intuitively reflects the clustering process. Coefficients of variation (CV) can reflect the degree of dispersion of data. It is generally considered that CV ≤ 0.1 shows a weak, CV > 0.1 ≤ 1 moderate, and CV ≥ 1 strong variation.

Factor Analysis (FA)
R-mode FA including Principal Components Analysis (PCA), was carried out to reveal the internal structure of the correlation matrix between variables. In general, this attempts to simplify the complex and diverse relationships that exist among a set of variables (measured geochemical constituents) by revealing common and unobserved factors that link together with the seemingly unrelated variables [36,41,42]. The PCA technique extracts the eigenvalues and eigenvectors from the covariance matrix of original variables. The common factors are the uncorrelated (orthogonal) variables, obtained by multiplying the original correlated variables with the eigenvector, which is a list of coefficients (loadings or weightings) [31,35,43]. The number of factors retained based on the Kaiser criterion, from which only the factors with eigenvalues greater than 1 are extracted [44]. The loading matrix of the original variables and the common factors are obtained by the method of varimax rotation, in which each factor is described only in terms of those variables and affords greater ease for interpretation [45]. The largest loading, either positive or negative, suggests the meaning of the dimension; positive loading indicates that the contribution of variables towards the geochemical processes increases with the increasing loadings in a dimension; whereas negative loading indicates a decrease [33,46,47]. Cluster analysis and factor analysis are performed with the IBM SPSS 22 statistical software.

Geochemical Modelling (GM)
Geochemical modelling is a powerful technique for characterizing geochemical phenomena and predicting their evolution in time, as well as in space when coupled with widely applied flow modelling [48][49][50]. PHREEQC, designed by the USGS, is one of the most popular software packages for modelling [51]. It is used here to calculate the saturation index of minerals and to conduct reverse modelling for the water-rock interactions on flow paths. The saturation index can qualitatively reflect the dissolution and precipitation of minerals during geochemical evolution. Based on the principle of mass conservation, inverse modelling calculates the molar amounts of minerals and gases that dissolve or precipitate (degas), and quantitatively explains the difference of hydrochemical components between initial solution and final solution, assessing the evolution of groundwater characteristics along the flow path [25,32]. The reaction path is determined by the hydrogeological conditions (hydrochemical characteristics, including changes in regional head and the geographical locations of sampling points) [34,52]. The determination of hydrogeochemical reactions along the flow path mainly depends on the determination of "possible mineral phases", and is mainly based on the lithologies and mineral compositions of aquifer host media, and on the chemical compositions and storage conditions of the groundwater [21].

General Hydrochemical Characteristics
Sampling information, chemical analyses, and hydrochemical types of the hot waters in the study area are shown in Table S1. Descriptive statistical analysis of the major and minor hydrochemical parameters of the 28 hot water samples is listed in Table S2. The pH of the water ranges from 6.2 to 8. This indicates that the hydrochemical characteristics of the hot water are controlled by rock weathering,

Hierarchical Cluster Analysis (HCA)
Based on the 11 main hydrochemical parameters (T, pH, K + , Na + , Ca 2+ , Mg 2+ , Cl − , SO4 2− , HCO3 − , SiO2, and TDS), Q-mode HCA was carried out on results from the thermal groundwater. The results of cluster analysis are represented using a dendrogram ( Figure 6). Classification is clearest when the distance is 10 ( Figure 6) and the 28 samples can be divided into three clusters: A, B, and C.

Hierarchical Cluster Analysis (HCA)
Based on the 11 main hydrochemical parameters (T, pH, K + , Na + , Ca 2+ , Mg 2+ , Cl − , SO4 2− , HCO3 − , SiO2, and TDS), Q-mode HCA was carried out on results from the thermal groundwater. The results of cluster analysis are represented using a dendrogram ( Figure 6). Classification is clearest when the distance is 10 ( Figure 6) and the 28 samples can be divided into three clusters: A, B, and C.  distance is 10 ( Figure 6) and the 28 samples can be divided into three clusters: A, B, and C. Observations of the dendrogram, provides information on the degree of similarity of the three categories. The link distances between clusters A and B are relatively small (16), indicating that the hydrogeochemical characteristics of these clusters are similar, whereas the link distances between cluster C and the others are up to 25. Thus, the hydrogeochemical characteristics are quite different from those of the other clusters. The Piper diagram (Figure 7) shows that the waters in all three clusters are of SO 4  Observations of the dendrogram, provides information on the degree of similarity of the three categories. The link distances between clusters A and B are relatively small (16), indicating that the hydrogeochemical characteristics of these clusters are similar, whereas the link distances between cluster C and the others are up to 25. Thus, the hydrogeochemical characteristics are quite different from those of the other clusters. The Piper diagram (Figure 7) shows that the waters in all three clusters are of SO4-Ca type. Among the clusters, the milliequivalent of SO4    (Table 1, Figure 8a,c,d,f,h). Sample H10-1 was the lowest Mg 2+ concentration and the highest HCO3 − concentration-80.4 mg/L and 195.9 mg/L-respectively. Based on the stratigraphy lithology of the study area, it can be inferred that cluster A waters may dissolve more evaporite (gypsum and anhydrite) and silicate rocks, but relatively less carbonate. Cluster B contains 14 samples, accounting for 50% of the total. The concentrations of each parameter in this cluster lies between those of the other two, with the exception of Na + and Cl − , where the concentrations of these ions are the maximum of all the water samples, at 37.0 mg/L and 36.0 mg/L, respectively (Table 1, Figure 8b,g). This may reflect the dissolution of more halite. Cluster C contains only two samples, accounting for about 7.1% of the total. However, the concentration of HCO3 − and pH of these samples are the highest (221.6 mg/L and 7.6, respectively), of all samples, whereas other parameters are the smallest, indicating that more carbonate was dissolved (Figure 8e,i). Thus, cluster A has high TDS, high Ca 2+ and SO4 2− , cluster B is medium TDS, and high Na + and Cl − , and cluster C has low TDS and high HCO3 − .  Figure 7. Piper diagram of the thermal groundwater samples.  Comparison of the TDS, Ca 2+ , SO4 2− , and HCO3 − contents in the three clusters of water samples show that TDS, Ca 2+ , and SO4 2− gradually decrease, whereas HCO3 − gradually increases. According to the dissolution order of anions in crustal rocks (Cl − > SO4 2− > HCO3 − ), and combined with the hydrochemical characteristics of the salty hot spring (Cl-Na type, TDS of 13.37 g/L) in Wenquanzhenn, eastern Sichuan basin, which is formed in the same geothermal reservoir [3], the geothermal groundwater in the study area (SO4-Ca type) is undergoing gradual desalination, but it has not yet reached the stage of 'fresh water' dominated by HCO3 − .

Factor Analysis
The correlation coefficient matrix is the basis of factor analysis. Table 2 illustrates the statistics of correlation among the hydrochemical parameters of hot water in the study area. There is a significant correlation among SO4 2− , Cl − , Ca 2+ , Mg 2+ , Na + and TDS (correlation coefficients are greater than 0.7), which indicates that there is information overlap among variables, and data dimensionality reduction analysis is needed.

Factor Analysis
The correlation coefficient matrix is the basis of factor analysis. Table 2 illustrates the statistics of correlation among the hydrochemical parameters of hot water in the study area. There is a significant correlation among SO 4 2− , Cl − , Ca 2+ , Mg 2+ , Na + and TDS (correlation coefficients are greater than 0.7), which indicates that there is information overlap among variables, and data dimensionality reduction analysis is needed. PCA is used to extract the eigenvalues from physico-chemical parameters of the 28 water samples in the study area. The eigenvalues of the extracted factors, their percentage of variance, and cumulative percentage of variance of the hydrochemical parameters are presented in Table 3. In particular, four factors are extracted based on the Kaiser criterion, which explain 83.86% of the total variance. These are considered to reflect the basic information of the original data. The loadings for varimax rotated factor matrix in the four factors model are compiled in Table 4. It is considered that the corresponding hydrochemical components are closely related to the hydrogeochemical processes represented by the factor when the absolute value of factor loading is more than 0.5 (marked in bold). The scores of each factor reflecting the degree of influence on each water sample are calculated, and the scatter plots of factor scores are drawn (Figure 9).  variance. These are considered to reflect the basic information of the original data. The loadings for varimax rotated factor matrix in the four factors model are compiled in Table 4. It is considered that the corresponding hydrochemical components are closely related to the hydrogeochemical processes represented by the factor when the absolute value of factor loading is more than 0.5 (marked in bold). The scores of each factor reflecting the degree of influence on each water sample are calculated, and the scatter plots of factor scores are drawn (Figure 9). Extraction method: Principal component analysis.  Table 4). The high value points of factor 1 scores are mostly from cluster A, indicating that in this cluster is greatly affected by Factor 1 (Figure 9a). The scores of samples corresponding to Factor 1 gradually decrease in the three clusters, implying that the effects of Factor 1 are gradually weakened. There is a good linear correlation between SO 4 2 and Ca 2+ in the geothermal water (R 2 = 0.92, Figure 10a), indicating that the dissolution of gypsum and anhydrite in the T 1 j 2 and T 2 l rocks is the main source of these components. Figure 10a shows  (Figure 10b,c), suggesting that the dissolution of carbonate rocks is not the main source of Ca 2+ and Mg 2+ . The dissolution of gypsum also increases the content of Ca 2+ in the hot water and inhibits the dissolution of carbonate rocks due to the common ion effect, consistent with the negative factor loading of HCO 3 − (−0.57). Mg 2+ may originate from the dissolution of magnesium sulfate minerals in the geothermal reservoirs, indicated by the relatively good correlation between Mg 2+ and SO 4 2 (R 2 = 0.69) ( Figure 10d). Factor 1 therefore, reflects the dissolution of gypsum and anhydrite, magnesium sulfate minerals and pyrite-bearing formation rocks. suggesting that the dissolution of dolomite, calcite, gypsum, and anhydrite is responsible for the production of Ca 2+ , Mg 2+ , SO4 2− and HCO3 − , thus eliminating the cation exchange between Na + , Ca 2+ , and Mg 2+ [50,54,55]. The incongruent dissolution of albite may contribute to the excess Na + . Thus, Factor 3 reflects the dissolution of halite and albite. Factor 4 corresponds to a variance portion of 10.8% with positive loadings for pH (0.62) and K + (0.76). The high value points of the scores are mainly from clusters A and B, and the trend is similar to that of Factor 1. The Factor 4 score gradually decreases from cluster A to cluster C, indicating that the influence of this factor gradually weakens. The lithology and factor analyses results indicate that the main source of K + is in dissolution of silicate rocks (K-feldspar). Thus, Factor 4 therefore reflects the dissolution of K-feldspar related to pH.

Saturation Index
The saturation index is an important measure, evaluating the equilibrium between water and minerals. Characteristics variation can be used to differentiate the stages of hydrochemical evolution and identify the main water-rock interactions controlling the process. The saturation indices of the 28 water samples and minerals in the study were calculated using PHREEQC software and are presented in Table S3. The results show that in all samples the saturation index of quartz is greater than 0 (0.04-0.62), with an average of 0.35, and that for chalcedony is also close to 0, with an average value of −0.024 (−0.32-0.24). Compared to quartz, chalcedony is closest to the saturation state, indicating that chalcedony most likely controls the SiO2 content of the water. Except for W16-1 and W9-1, the saturation index of anhydrite for the water samples is less than 0 (Figure 11a), and thus gypsum is undersaturated for most of the samples (Figure 11b). Anhydrite and gypsum are wellcorrelated with (Ca 2+ + SO4 2− ) (R 2 is 0.70 and 0.99, respectively), indicating that the dissolution of suplhate minerals is the main source of Ca 2+ and SO4 2− . Calcite and dolomite are both supersaturated in most of the water samples (Figure 11c-e), confirming that the dissolution of carbonates is limited by the increase in the Ca 2+ content of the water due to the dissolution of gypsum. Halite is  (Table 2), and SiO 2 geothermometers are widely used to estimate the temperature of reservoirs waters. The positive loadings for HCO 3 − and SiO 2 indicate that the increasing SiO 2 content of geothermal waters may promote an increase in HCO 3 − concentration. Factor 2 reflects the dissolution of quartz or chalcedony and other silicates. Factor 3, with high positive loadings for Cl − (0.89) and Na + (0.93), accounts for 16.76% of the total variance. The high value points of the scores are mainly from cluster B, which shows that these water samples are greatly affected this factor (Figure 9b). The ionic ratios of Cl − and Na + (Figure 10e) show that most of the water falls above the dissolution line for halite, indicating that Na + has other sources apart from the dissolution of halite. The ratio of (Ca 2+ + Mg 2+ )/(SO 4 2− + HCO 3 − ) is close to 1 (Figure 10f), suggesting that the dissolution of dolomite, calcite, gypsum, and anhydrite is responsible for the production of Ca 2+ , Mg 2+ , SO 4 2− and HCO 3 − , thus eliminating the cation exchange between Na + , Ca 2+ , and Mg 2+ [50,54,55]. The incongruent dissolution of albite may contribute to the excess Na + . Thus, Factor 3 reflects the dissolution of halite and albite. Factor 4 corresponds to a variance portion of 10.8% with positive loadings for pH (0.62) and K + (0.76). The high value points of the scores are mainly from clusters A and B, and the trend is similar to that of Factor 1. The Factor 4 score gradually decreases from cluster A to cluster C, indicating that the influence of this factor gradually weakens. The lithology and factor analyses results indicate that the main source of K + is in dissolution of silicate rocks (K-feldspar). Thus, Factor 4 therefore reflects the dissolution of K-feldspar related to pH.

Saturation Index
The saturation index is an important measure, evaluating the equilibrium between water and minerals. Characteristics variation can be used to differentiate the stages of hydrochemical evolution and identify the main water-rock interactions controlling the process. The saturation indices of the 28 water samples and minerals in the study were calculated using PHREEQC software and are presented in Table S3. The results show that in all samples the saturation index of quartz is greater than 0 (0.04-0.62), with an average of 0.35, and that for chalcedony is also close to 0, with an average value of −0.024 (−0.32-0.24). Compared to quartz, chalcedony is closest to the saturation state, indicating that chalcedony most likely controls the SiO 2 content of the water. Except for W16-1 and W9-1, the saturation index of anhydrite for the water samples is less than 0 (Figure 11a), and thus gypsum is undersaturated for most of the samples (Figure 11b

Reverse Modelling of Hydrochemical Pathways
Determination of the modelling path: The genetic model of the thermal groundwater in the Chongqing area is summarized as follows [4,29]: Groundwater receives recharge from precipitation in areas of outcrop of carbonate rocks in the cores of the anticlines. The cool water flows through these rocks along the flanks of the anticlines to deeper areas. After being heated by the geothermal heat flow, the hot water rises and issues in the form of hot springs at low points in the topography at the middle or the plunging ends of the anticlines, appearing in wells in the flanks of the structures. Rainwater (RW) was selected as the initial water for the reverse modelling in this study, with composition taken from the average of nine rainfall events from 2009 to 2011 in the Beibei District of

Reverse Modelling of Hydrochemical Pathways
Determination of the modelling path: The genetic model of the thermal groundwater in the Chongqing area is summarized as follows [4,29]: Groundwater receives recharge from precipitation in areas of outcrop of carbonate rocks in the cores of the anticlines. The cool water flows through these rocks along the flanks of the anticlines to deeper areas. After being heated by the geothermal heat flow, the hot water rises and issues in the form of hot springs at low points in the topography at the middle or the plunging ends of the anticlines, appearing in wells in the flanks of the structures. Rainwater (RW) was selected as the initial water for the reverse modelling in this study, with composition taken from the average of nine rainfall events from 2009 to 2011 in the Beibei District of northwest Chongqing [56]. The average values of the three groups obtained by cluster analysis are used as the final water. Three modelling paths are selected: defining the evolution of rainwater to Clusters A, B, and C.
Determination of possible mineral phases: The dissolution and precipitation of calcite, dolomite and gypsum in the geothermal reservoir controls the proportions of Ca 2+ , Mg 2+ , SO 4 2− and HCO 3 − in the geothermal waters. Changes in Cl − and Na + concentrations are attributed to the dissolution of halite, and the incongruent dissolution of albite may be another source of Na + . The dissolution of K-feldspar and chalcedony may explain changes in K + and SiO 2 content, respectively. Changes in dissolution and evasion of CO 2 may also occur as the waters evolve. Calcite, dolomite, gypsum, halite, albite, K-feldspar, chalcedony, kaolinite, CO 2 , and water were selected as reaction phases, and the mineral transfers obtained by inverse modelling using PHREEQC are shown in Table 5. Discussion of modelling results: Results show that similar hydrochemical reactions took place in the three reaction paths. Albite, dolomite, gypsum, halite, K-feldspar, and CO 2 all dissolve and enter the hot water whereas calcite, chalcedony and kaolinite all precipitate. The implied dissolution reactions are shown in Table 6 The dissolution of halite is 0.3196 mmol/L, increasing the concentrations of Na + and Cl − in the water, and the dissolution of albite (0.5746 mmol/L) is an important additional source of Na + . The change in K + concentration is mainly attributed to the dissolution of K-feldspar (0.5719 mmol/L). Model 1 shows that kaolinite precipitates (−0.5733 mmol/L), mainly as a result of dissolution of albite and K-feldspar (Table 5). Modelling results also show the precipitation of chalcedony, contrary to the results of the dissolution of chalcedony noted above. This is because only the mass balance of Si is taken into account in the modelling and the precipitation condition of SiO 2 (reaching the solubility of amorphous SiO 2 ) is neglected. The dissolution of both albite and K-feldspar increase the content of HCO 3 − in the groundwater, consistent with the results of factor analysis. Comparison of the amount of minerals dissolved in the three reaction pathways indicates that the dissolution of gypsum gradually decreases, from Model 1 to Model 3 (19.57 mmol/L, 17.18 mmol/L and 10.42 mmol/L), consistent with the decreasing trend in concentrations of Ca 2+ and SiO 4 2− from cluster A to cluster C. The dissolution of dolomite also decreases, resulting in a decrease in HCO 3 − in the water. The dissolution of CO 2 and silicate minerals and the precipitation of calcite leads to an increase in HCO 3 − , producing a gradual increase in concentrations in the three clusters. The dissolution of K-feldspar decreases from Model 1 to Model 2, whereas halite and albite show increasing trends with a significantly increasing dissolution of halite resulting in the concentration of Na + in cluster B, significantly greater than that in cluster A. However, the dissolution of halite and albite gradually decreases from Model 2 to Model 3. In summary, the dissolution of mineral phases in the Chongqing geothermal groundwater gradually decreases as it evolves. Compared to the hydrochemical characteristics of the salty hot spring (Cl-Na type, TDS of 13.37 g/L) in Wenquanzhenn, eastern Sichuan basin, which is formed in the same geothermal reservoir [3], the hot water is in the stage of desalination, but has not reached the stage of fresh water dominated by HCO 3 − .

Conclusions
The distribution and characteristics of the thermal groundwaters in Chongqing are strongly controlled by the anticlinal structures in the area. The pH of the water ranges from 6.2 to 8.5, and this low-to-middle temperature hot water (32.9-57 • C) shows a slight cooling with time (2009-2016) most likely influenced by their degree of development and utilization. The groundwater is of SO 4 -Ca type with the TDS ranging from 1620 to 2929 mg/L. It can be conjectured that for highly exploited hot waters, TDS gradually decreases; whereas less exploited waters show an increase. The major ions are Ca 2+ , Mg 2+ , SO 4 2− , and HCO 3 − , but waters are poor in Na + , K + , and Cl − , with only trace amounts of SiO 2 . Both Li + and Sr 2+ are relatively abundant in the hot water. The three clusters of water samples are summarized using hierarchical clustering analysis (HCA): high TDS-high Ca 2+ and SO 4 2− water (cluster A); medium TDS-high Na + and Cl − water (cluster B); low TDS-high HCO 3 − water (cluster C), respectively. Combined with previous analyses, data provide a clear understanding of the origin and evolution of major ions in the thermal groundwater. Ca 2+ and SO 4 2− mainly derive from the dissolution of gypsum and anhydrite in the Triassic rocks but, in addition, SO 4 2− also increases by dissolution of pyrite in coal-bearing deposits. Mg 2+ may be derived from the dissolution of magnesium sulfate minerals (epsomite) in the reservoir, resulting in the relatively good correlation between Mg 2+ and SO 4 2− (R 2 = 0.69). Although calcite is precipitated, the dissolution of dolomite and CO 2 are the main sources of HCO 3 − , and the dissolution of dolomite also releases Mg 2+ into the hot water. The dissolution of halite is the main cause of the increase in Na + and Cl − . The dissolution of CO 2 promotes the dissolution of albite and K-feldspar, adding Na + and K + to the water. The saturation index shows that the characteristic component (SiO 2 ) of the water is controlled by the dissolution of chalcedony. Taking the rainwater of the northwest of the study area as the initial water and three clusters of water samples as the final water the hydrochemical pathway modelling, quantitatively analyzes the water-rock interaction during the evolution of the water. Comparing the mineral transfers in the three models, it is found that-except for halite and albite in Model 2-the dissolution of all minerals decreases. Combined with the hydrochemical characteristics of hot water in the same reservoir in the adjacent area (Cl-Na type, TDS of 13.37 g/L), we believe that the thermal groundwater in Chongqing is in the stage of desalination, but it has not yet reached the stage of fresh water dominated by HCO 3 − .
In this paper, traditional graphic methods, multivariate statistics and hydrogeochemical modelling are used to discuss the hydrochemical characteristics, source and evolution of dissolved components in the thermal groundwaters of the Triassic carbonates in Chongqing, China. Distribution of the waters is controlled by anticlines and the study provides a deeper and more comprehensive knowledge for the genesis mechanism and evolution trend of SO 4 -Ca type water in carbonate reservoirs in an exploiting environment.