Promissory geothermal zones of the Sonora State

The study area containing the promissory geothermal zones is located approximately 50 km in straight line to the east of Hermosillo city (the capital of the Sonora State), in the SMO mountain portion of east-central Sonora (see Figures 1b and 1c). Schematic geological cross sections of the study area are shown in Figure 3. In this area, twelve geothermal localities were identified with a hydrothermal activity in association with hot spring waters (see Figure 1c). These localities have been geochemically grouped in four major hydrothermal zones: the Northwest (NW) zone which is characterized by the shallow hot springs from Aconchi (ACH) and Cumpas (CMP); the Northeast (NE) zone which is characterized by the shallow hot springs from Tonibabi (TNB), Huasabas (HSB), Divisaderos (DVS), Granados (GRN), and Bacadehuachi (BCD); the Central (C) zone which is characterized by the shallow hot springs from Matape (MTP) and Arivechi (ARV); and the South zone (S) which is characterized by the shallow hot springs from San Marcial (SMR), Tecoripa (TCP) and Tonichi (TCH).

Figure 3 Schematic geological cross sections of the northeastern and southern zones of the study area (same legend as in Figure 2). 

The NE zone contains a Lower Cretaceous detritus-carbonated sedimentary sequence (^{Lawton et al., 2004}), which is cut by Upper Cretaceous volcanic rocks, and correlated with the Tarahumara Formation. Both sedimentary and volcanic sequences are intruded by a series of plutonic rocks (^{Roldán-Quintana, 1994}; ^{Almirudis, 2010}) (Figure 2). Tertiary rocks correlated with the Upper Volcanic Complex cover the whole series of rocks (^{Cochemé, 1985}). The basins are filled with alluvium and sediments from the Baucarit Formation, and, particularly the Moctezuma valley, with basalts of the Quaternary Moctezuma volcanic field. The NW, C, and S zones present similar lithologies to those in the NE zone, although lacking Quaternary volcanism. Furthermore, some of their outcrops correspond to ancient rocks composed mainly by carbonated sedimentary rocks (limestones and dolomites) and sandstones (see Figures 2 and 3).

Climatological data

The climate of the study area can be categorized as semiarid. The highest annual rainfall generally occurs on July, and it ranges from 114 to 175 mm, whereas for the dry season, which generally occurs on May, the rainfall is lower than 3 mm. The lowest and highest annual evaporation values are 84 mm and 331 mm, respectively. The air temperature over this area varies from a maximum value of 42 °C in the summer time (on June) up to a minimum temperature of 2 °C in winter (mainly on January). Table 1 summarizes the average annual climatological measurements collected from six meteorological stations located in the study area (BCD, HSB, DVS, MTP, ARV, and TCH). These data were compiled from the Eric III Climatological Database v.1.1 (^{IMTA, 2013}).

Fluid sampling and chemical analyses

For the fieldwork, a geochemometric programme was used for the collection and analysis of representative water samples from the identified hydrothermal zones. Major physicochemical parameters were in situ measured. A HACH portable multimeter (equipped with temperature and electrical conductivity sensors), and a HANNA portable potentiometer were used for measuring temperature, and pH of the hot spring waters. All the pH measurements were recorded with the temperature compensation option after calibration with standardized buffer solutions of pH 4.0, 7.0 and 10.0. Measurement errors of pH (accuracy and precision) were typically better than ±0.1 pH units (±1.4%); whereas for the temperature and electrical conductivity measurements, the errors were better than ±1.1 °C (±2.6%) and ± 33 µS/cm (±3.0%), respectively. Ten replicates of all these measurements were systematically recorded in the hot spring discharges. Total alkalinity (HCO₃ ^- and CO₃ ^2-) was determined in situ, by triplicate, after sampling by using a titration method with a micro-burette according to the standardized analytical methodology proposed by ^{Nicholson (1993)}. Gas samples were not collected because steam vents or fumaroles were not observed in the study area. Small gas bubbling sites were intermittently found in a few zones with very low discharge flow rates, which made difficult their collection.

Polyethylene bottles of 500, 100, and 50 mL volume (previously washed with HCl and rinsed three times with deionized water) were used for collecting water samples to be used in the chemical analysis of major (anions and cations, and other compounds) and trace (metals) elements, as well as for the isotopic analyses (¹⁸O/¹⁶O, and D/H). For the analysis of anions, hot spring waters were collected without adding acid, whereas for the determination of cations, the water samples were preserved by adding ultrapure HNO₃ to pH < 2. All these water samples were previously filtered through a 0.2 μm membrane filter.

Chemical analyses (major and trace elements) were carried out at Actlabs Laboratories (Canada), whereas the isotopic analyses were conducted at the Isotope Laboratory of the University of Arizona (USA). Major cations and trace elements were determined by Inductively Coupled Plasma-Mass Spectrometry and Inductively Coupled Plasma-Optical Emission Spectrometry, whereas the anions were analysed by ion chromatography. The stable isotopes of ¹⁸O/¹⁶O and D/H were analysed by gas source Isotope-Ratio Mass Spectrometry. Table 2 summarizes the precision and accuracy errors (in %) reported for the chemical analyses, including the limits of detection (LOD’s), for each one of the elements. Accuracy was evaluated by analysing water certified standards as external testers. Ionic charge balances were calculated among the major cations (Li⁺, Na⁺, K⁺, Ca²⁺, and Mg²⁺) and anions (F⁻, Cl⁻, NO₃ ^-, SO₄ ^2-, PO₄ ^3-, and HCO₃ ^-) for checking the quality of the chemical analyses of hot spring waters.

Hydrogeochemistry

Solute geothermometry and mixing models

Reservoir temperatures were estimated by solute geothermometry using the SolGeo software, which was programmed by ^{Verma et al. (2008)}. The following solute geothermometers were used in this study: (i) Na-K geothermometers: TNaKF79 (^{Fournier, 1979}); TNaKVS97 (^{Verma and Santoyo, 1997}); and TNaKDSR08 (^{Díaz-González et al., 2008}); (ii) Na-K-Ca geothermometer: TNaKCaFT73 (^{Fournier and Truesdell, 1973}); (iii) Na-Li geothermometer: TNaLiVS97 (^{Verma and Santoyo, 1997}); and (iv) SiO₂ geothermometer: TSiO2VS97 (^{Verma and Santoyo, 1997}).

Table 3 summarizes the geothermometer equations and their applicability constraints. The Na-K geothermometers were selected for the prediction of the equilibrium temperatures because these tools generally provide more consistent and reliable temperature estimates than other solute geothermometers (^{Verma et al., 2008}). The K-Mg, Na-K-Ca, Na-Li, and SiO₂ geothermometric equations have been used for obtaining low and high temperatures of the system, keeping in mind that these tools should be applied with caution.

To evaluate the actual equilibrium state that exhibit the fluids collected in the hydrothermal zones of Sonora, the well-known Na-K-Mg ternary diagram proposed by ^{Giggenbach (1988)} was used. The equilibrium isothermal curves were plotted from 75 °C to 350 °C by using the equations of the Na-K and K-Mg geothermometers (TNaKVS97 and TKMgG88, respectively). The hot spring waters sampled in the geothermal system were also plotted in the same ternary diagram.

For ascending hot spring waters that lack equilibrium conditions (i.e., temperatures and chemical compositions falling out of the full equilibrium curve of the Giggenbach’s ternary diagram), the influence of conductive cooling processes was examined by using a well-known mixing geochemical model (i.e., a plot between the silica and fluid enthalpy) to evaluate the mixing degree with colder groundwaters assuming slow fluid flow rates (usually <0.5 L/s) and lengthy flow paths (^{Nicholson, 1993}). ^{Fournier (1977)} stated that, whether or not equilibrium is achieved after mixing, the equilibrium temperature of hot spring waters may not be directly estimated from either a simple solubility relationship or solute geothermometer unless the mixing process is analysed. The application of mixing models to hot spring and cold waters by using silica composition and temperatures was performed for inferring the original fluid composition and temperature at reservoir conditions. In these models, the initial silica composition of hot spring waters is inferred from the curve of quartz solubility (with an approximated uncertainty of ±5%), assuming that no further solution or deposition of silica occurs before or after mixing.

Fluid-mineral equilibrium

For this geochemical study, it was important to understand how the hot spring waters chemically interact with host rock-minerals (i.e., fluid fractionation and chemical mobility processes), as well as their transient transport towards the surface. The local equilibrium state between the hydrothermal waters and the surrounding mineral phases was thermodynamically investigated by using a comparison between the water composition and the solubility of specific minerals at the most probable reservoir temperature conditions (^{Torres-Alvarado, 2002}). The prediction of the distribution of chemical species in the fluid, and the saturation indexes of minerals were then determined by using the software “The Geochemist’s Workbench” (GWB v. 8.0: React and Act2 programs). React is a computer program for modeling equilibrium states and geochemical processes in systems that contain an aqueous fluid; whereas Act2 is a program that calculates and plots activity-activity diagrams showing the stability of minerals and the predominance of aqueous species in geochemical systems.

React program calculates the equilibrium distribution of aqueous species in a fluid, a fluid’s saturation state with respect to minerals. The program may also trace the geochemical evolution of a system as it undergoes reversible or irreversible reaction in an open or closed system, either at a given temperature or for an interval of temperatures. GWB software uses a thermodynamic database created by the Lawrence Livermore National Laboratory (^{Delany and Lundeen, 1990}). As part of this modelling, a set of non-linear governing equations describing the equilibrium state of a geochemical system is generally proposed. The resulting non-linear equation system is solved for multi-component equilibrium by means of the Newton-Raphson method. A full description of the mathematical formulation is reported by ^{Bethke (2008)}.

For plotting the ionic activity diagrams of the hot springs, the following assumptions were considered: (i) an aqueous solution should be continuously present in the system; (ii) the aluminium (Al³⁺) content in the minerals must be present as a solid phase; (iii) the silica concentration must be determined by the solubility curve of quartz; and (iv) the temperature and the corresponding pressure of the system are constant. The ionic activities in solution corresponding to the water samples at pseudo- or equilibrium conditions were obtained as activity-activity diagrams. The results inferred from this fluid-equilibria modelling were examined with caution because the steam/gas phase composition for the fluid samples was not available. The lack of gas-phase composition limited the use of the multi-component geothermometry methodology [i.e., log(Q/K) curves depending on temperature] for inferring the equilibrium temperatures with accuracy. Regarding this missing information, ^{Reed and Spycher (1984)} and ^{Pang and Reed (1998)} identified several uncertainties that may affect the calculation of equilibrium temperatures, among which are the following sources: (i) the large pH differences expected between the fluid with and without the gas-phase composition, which generally result in large uncertainties for the calculation of the total molalities and ion activities; (ii) the effect of degassing on the aqueous species, which in combination with pH have significant effects for mineral super- or under-saturation; (iii) the effect of boiling may also produce an irregular dispersion of log(Q/K) curves due to the combined effects of aqueous component concentration caused by water loss, and the pH change due to CO₂ loss; (iv) the effect of dilution also produces a shift and dispersion of the mineral log(Q/K) curves causing a displacement in most curve intersections (with log Q/K = 0) with implications on the calculation of low equilibrium temperatures; and (v) the uncertainty of the aluminium analyses (or analytical errors) which may increase the propagated errors by up to 35 ^oC in the estimation of reservoir equilibrium temperature.

Fluid sampling

From the fluid sampling campaign, it was observed that the ACH, TNB, TCP, and TCH hot spring waters are geologically correlated with the Laramide granitoids (Figure 3), and probably associated to NW-SE oriented regional faults at the basin limits. The CMP and DVS hot spring waters rise towards the surface probably through fractures and secondary faults and interact with ignimbrite outcrops characteristic of the Sierra Madre Occidental (Figure 3). The BCD, GRN, and SMR thermal spring waters ascend to the surface directly on alluvium, whereas the HSB hot springs ascend on ancient alluvium that has been correlated with the Baucarit Formation (Figure 3). The MTP hot spring water upflow to the surface occurs perhaps through a more than 10 km long regional fracture present in limestones and covered with sediments, whereas the waters of hot spring ARV rise towards the surface to interact with limestones from older Formations.

Physicochemical measurements

The physicochemical parameters measured in situ in the hot springs are summarized in Table 4. The highest shallow temperatures, which ranged from 41.5±7.5 °C to 63.1±0.6 °C, were measured for hot spring waters that emerge mostly on granite formations: TCP (S), GRN (NE), TNB (NE), ACH (NW), and TCH (S). Intermediate shallow temperatures between 31.4±0.2 °C and 40.0±0.2 °C were measured for thermal springs emerging on sedimentary formations: SMR (S), BCD (NE), HSB (NE), and ARV (C). The lowest temperatures were measured in waters coming from the hot springs DVS (29.3±0.2 °C), MTP (30.9±2.9 °C), and SMR (31.4±0.2 °C), which emerge on ignimbrite, limestone, and alluvium rocks, respectively.

In the case of the pH measurements, the values varied from 6.0±0.1 (for the DVS hot spring) to 9.0±0.1 (for the HSB), whereas those related with the electrical conductivity show a larger variability with values ranging from 426±12 μS/cm (for the BCD hot spring) to 3670±110 μS/cm (for the DVS). Total dissolved solids (TDS) were also estimated for all the hot springs sampled, and they ranged between 343±10 ppm (for the TCP hot spring waters, S) and 4519±136 ppm (for the DVS, NE): see Table 4.

Shallow temperature, pH, and electrical conductivity measurements were plotted in distribution maps using the Surfer software (Figure 5). These distribution maps show the spatial distribution of these parameters and the following geochemical signatures: (i) the hot spring waters with shallow temperatures greater than 40 °C are mostly found in the NW, NE and S hydrothermal zones, whereas the lower temperatures were mostly observed at the C zone, and the anomalous saline hot spring DVS (NE); (ii) a pH transition from 6.8 (neutral) to 9.0 (strongly alkaline) was observed for the fluids from the C to NE geothermal zone, with the exception of the anomalous DVS sample which exhibits a pH of 6.0 (slightly acidic); and (iii) the highest ionic concentrations (given by the high electrical conductivity values) were located at the NE zone, whereas the lowest values were mostly found at the C zone, which agree with the temperature and pH patterns observed probably as a product of the WRI processes.

Figure 5 Distribution of physicochemical parameters measured in the hot spring waters of the Sonora geothermal system: (a) temperature, (b) pH, and (c) electrical conductivity. 

Chemical analyses of water samples

The chemical analyses for the hot spring waters sampled in the geothermal system of Sonora under study are compiled in Table 5. The content of major elements (anions and cations) and SiO₂ (in mg/kg) and the isotopic analyses of δ¹⁸O and δD (in ‰) have been reported in the first section of Table 5. The next section was used to report the abundance of trace elements in µg/kg, whereas the last section was used to present the concentration of the rare earth elements (REE) in µg/kg. A statistical descriptive plot showing the main concentration patterns of major and trace elements in all hot spring samples is presented in Figure 6 (as box-plots with some statistical parameters given by the median, percentile, and min-max values), whereas the REE average concentration patterns for the four hydrothermal zones (normalized with respect to chondrite) are plotted in Figure 7.

Figure 6 Statistical descriptive plot showing the main compositional features observed for major and trace elements measured in the hot spring samples of the Sonora geothermal system. 

Figure 7 Rare Earth Element (REE) patterns normalized to chondrite values for the hot spring waters of the NE, NW, C and S zones. All the symbols correspond to the mean value for each zone each zone; bars indicate minimum and maximum values. 

The statistical descriptive plot (Schoeller type) clearly indicates that the hot spring waters are mainly depleted in the major elements F^-, Cl^-, K⁺ and Li⁺, and enriched in HCO₃ ^-, SO₄ ^2-, Na⁺ and Ca²⁺ (Figure 6). Regarding the trace element contents, depletions were mainly observed for Cs, Sr, Sb, V, Co and Tl, whereas enrichments were mostly given for Rb, Ba, As, Al, Mn and Fe. A high chemical mobility from the host rock to the fluid (fractionation) was mainly observed in Mn, Fe, Sr, As, Rb, Ba, including the light REE (La and Ce). Intermediate and heavy REE were immobile elements that were not partitioned into the fluid due to the alkaline pH conditions (pH>6) that dominated most of the WRI and water mixing processes observed in the hydrothermal system.

An interesting geochemical signature of the NE hot spring waters is the zig-zag pattern observed for the REE, with positive and negative anomalies of Eu and Gd, respectively. This REE pattern is roughly observed for the NW and S hot spring waters, whereas for the C zone appears almost as a flat pattern (see Figure 7).

As seen in Table 5, boron (a highly mobile and conservative element) was not measured in the hot spring waters because it is generally present in high-temperature geothermal systems (^{Aggarwal et al. 2000}; ^{Gurav et al., 2016}). According to the local WRI processes (i.e., fluid fractionation or chemical mobility) that occur in the hydrothermal system of Sonora, the non-volatile elements characterized as non-reactive or conservative compounds were mainly associated to the contents of Cl, Li, Rb, and Cs, whereas the reactive or geoindicator elements were typified by the concentrations of major cations (Na, K, Ca, Mg), and SiO₂.

A multivariate statistical analysis of the linear relationships for all the physicochemical parameters, major and trace elements was applied by using a correlation matrix (Table 6a and 6b). Correlations among these variables provided valuable information of the major hydro-geochemical processes that controlled part of the chemical signatures found in the hot spring waters. Applying the methodology suggested by ^{Bevington and Robinson (2003)} for a small data sample (n=12), correlation coefficients of r>0.823 provide an acceptable correlation for a confidence level of 99.9%, whereas r values ranging from 0.5 to 0.7 shows a moderate correlation among the parameters. From this analysis, the physicochemical parameters (EC and TDS) resulted in high correlations (r>0.82) for both the cation (Na⁺, K⁺, Li⁺ and Ca²⁺) and anion (SO₄ ^2- and HCO₃ ^-) contents (see Table 6a). From the correlation matrix, it was found that TDS and EC of the hot spring waters were mainly controlled by Li⁺-HCO₃ ^- and Na⁺-SO₄ ^2- ions, respectively.

On the other hand, the cation concentrations (Na⁺, K⁺, Li⁺, and Ca²⁺) were highly correlated with most of the major elements found in the hot spring waters, as well as with the trace elements Mn, Fe, Ge, Rb, Y, Zr, and Nb. In relation to the anion compositions, SO₄ ^2- resulted in a good correlation with Be, Ti, Ge, and Rb; whereas the HCO₃ ^- concentrations were well correlated with Mn, Fe, Y, Zr, and Nb. Among trace elements, Mn, Fe, Y, and Nb showed the highest correlation coefficients.

Hydrogeochemistry (fluid classification)

Major components of the hot spring waters were plotted in a Piper plot (Figure 8). As seen from this diagram, the geothermal fluids produced in the hot springs of the Central zone (MTP and ARV) were classified as calcium and magnesium-bicarbonate (Ca-Mg-HCO₃) waters, compositions that are in concordance with the sedimentary formation from where they emerge after WRI and mixing processes. The waters produced in the NW (ACH and CMP), two NE (GRN and TNB), and S (TCP, SMR and TCH) hot springs were mostly grouped as sodium-sulphate (Na-SO₄) waters, whereas those fluids coming from three hot springs from the NE hydrothermal zones (BCD, HSB and DVS) were mainly characterized as sodium-bicarbonate (Na-HCO₃) waters. As seen in Figure 8, the hot spring waters of the Central zone (MTP and ARV) have entirely different chemistry compared to the other springs, which indicates that these waters probably come from a totally different fluid source or reservoir. The major-ion data of the remaining hot spring waters very well matches with more systematic chemical signatures of a geothermal system of low-to-medium temperature, which mostly dominates along the NW, NE and S zones.

Figure 8 Piper diagram for the hot spring waters showing their geochemical classification. 

An additional geochemical classification of fluids is also represented in the Cl-HCO₃-SO₄ ternary diagram proposed by ^{Giggenbach (1988)}. According to this diagram (Figure 9), the hot spring waters produced at the Central (MTP and ARV), South (TCP and SMR) and NE (BCD and DVS) zones may be roughly classified as peripheral waters (i.e., bicarbonate and diluted chloride waters with a low concentration of sulphates) with a tendency toward volcanic waters; whereas the remaining hot spring waters fall on the plot section of “steam heated waters” due to the relative high proportion of sulphates found in the samples (NW: ACH and CMP; NE: HSB, GRN and TNB; and the anomalous TCH sample of the S zone, which seems to be the most representative water sample by its higher surface temperature, and typified as a mix chloride-sulphate or volcanic water). Some water samples (e.g. GRN, TNB, and HSB) exhibited anomalous alkaline pH between 6.8 and 9, which could not be strictly considered as characteristic of “steam heated waters”. The classification of these samples is difficult to understand because it would be expected to have acid waters (pH<3) instead of slightly acid sulphate waters with pH >6. To support these anomalous geochemical signatures, some recent geochemical findings reported in the literature on the possible origin of these geothermal fluids are discussed.

Figure 9 Ternary diagram showing the variation of major anions (Cl^-, SO₄ ^2- and HCO₃ ^-) in the hot spring waters. 

^{Smith et al. (2010)} pointed out that in magmatic-hydrothermal and associated geothermal systems, acid magmatic fluids (with pH<3) are usually typified by a high concentration of sulphates, which are commonly discharged from hot springs nearby to the vent of active (degassing) volcanoes. However, these authors also describe that slightly acid geothermal fluids (with pH>5) may be found close to lateral outflows some distance away from the main degassing vent. Under these conditions, “alkaline sulphate waters” might be discharged with a diluted ionic composition.

Furthermore, reactions among dissolved carbon dioxide and host rocks tend to form HCO₃ ^-, the concentration of which is probably affected by permeability and lateral flow. As a result of this process, hot springs that are fed directly from the reservoir tend to exhibit the lowest HCO₃ ^- concentrations (which is actually characterized by the TCH water sample: 38.6 mg/kg and the highest surface temperature). This enables the HCO₃ ^-/SO₄ ^2- ratio to be used as an indicator of flow direction. The flow of a fluid away from the upflow yields the suitable conditions for rock-water reaction, and therefore for an increase of HCO₃ ^-. This, combined with the loss of H₂S by rock-water reactions with increased lateral flow, leads to an increase in the HCO₃ ^-/SO₄ ^2- ratios away from the upflow zone with values ranging from 0.07 (TCH) to 2 (TCP) for hot springs located at low altitudes of ~400 m a.s.l. (S zone), and from 0.25 to 11 for hot springs situated at an interval of altitudes between 600 and 900 m a.s.l. (NW, NE, and C zones).

Additionally, and to support the detection of relative high proportions of sulphates in six water samples from the NW, NE and S zones (see Figure 9), three precipitated salt samples were collected in situ at the CMP, TCH and TNB hot springs, and analysed by X-Ray Diffraction (XRD) for determining the major mineral compositions. The results of these XRD analyses show that the minerals that precipitated from the CMP waters were mainly characterized as quartz (SiO₂), thenardite (Na₂SO₄), and calcite (CaCO₃), whereas those for TCH and TNB hot springs were typified as gypsum (CaSO₄ ^.2H₂O), and a mixture of thenardite (Na₂SO₄) and halite (NaCl), respectively (see Figure 10). Taking into account these XRD analyses, the composition of these hot spring waters (as a product of the interaction with or dissolution of these precipitated minerals) may be in agreement with the unusual concentrations of sulphates reported for these water samples (see Figure 9). All these classification signatures are mostly in agreement with the grouping inferences obtained from the Piper diagram.

Figure 10 XRD diffractograms of precipitated salt samples collected at the (a) CMP, (b) TCH and (c) TNB hot springs. 

Ternary diagram of Li-Rb-Cs

The lithium concentrations in the hot spring waters varied between 0.66 and 3.1 mg/kg for the NE zone, whereas lower concentrations were found for the NW and S hydrothermal zones which ranged from 0.41 to 0.47 mg/kg, and from 0.17 to 0.53 mg/kg, respectively. Hot spring waters that emerge from the C zone show the lowest lithium concentrations (0.020 mg/kg). On the other hand, Rb concentrations were mostly the highest for the NE zone (up to 123 µg/kg), whereas for the NW and S zones, intermediate concentrations ranging from 25.40 to 35.80 µg/kg, and from 9.95 to 75.3 µg/kg, respectively, were found. The lowest concentrations of Rb were actually measured for the C zone. In relation to Cs, a systematic composition pattern was identified for all the hydrothermal zones (ranging from 4 to 123 µg/kg for the NE; from 12.50 to 28.80 µg/kg for the NW; from 3.16 to 16.4 µg/kg for the S; and from 2.67 to 3.16 for the C). This pattern is shown in the ternary diagram of Li-Rb-Cs (Figure 11). The plot shows chemical signatures related to a group of hot spring samples falling between the basalt and rhyolite rock boundaries, which are well characterized by K/Rb and K/Cs ratios ranging from 109 to 192, and from 112 to 567, respectively. According to ^{Goguel (1983)}, it is assumed that these waters may interact with illite at temperatures between 190 and 210 ^oC.

Figure 11 Ternary diagram showing relative Li, Rb, and Cs contents in the hot spring waters. Outflow and upflow are schematically shown. Dotted area comprises the group of hot spring samples falling between the basalt and rhyolite rock boundaries, which limit the waters thay may have interacted with illite at temperatures between 190 and 210 °C (^{Goguel, 1983}). 

On the other hand, the enrichment of Rb in presence of low concentrations of Li seems to indicate that these hot waters reached temperatures around 210 °C at the reservoir conditions before mixing with shallow cold waters during their ascent towards the surface. As a result of these chemical signatures, the study of the ternary diagram of Li-Rb-Cs may also be used as a secondary proxy indicator for inferring the theoretical reservoir temperatures of the Sonora hot spring waters.

Isotopic analyses

The water isotopic composition (D/H and ¹⁸O/¹⁶O) is considered as a conservative signature, and it is a good geochemical indicator for inferring water origin, mixing and evaporation processes, as well as a proxy for identifying the intensity degree of water rock interaction when the rock permeability is low (^{Bahati et al., 2005}). Results of the isotopic analyses (δ¹⁸O and δD) for the hot spring samples collected in Sonora are reported in Table 5. The isotopic composition of these waters shows an approximated dispersion of 2.7 units (from -6.9 to -9.6 with an accuracy of ±0.08 ‰) for ¹⁸O/¹⁶O, and 19 units (from -50 to -69 with an accuracy of 0.9 ‰) for D/H.

All these data points are shown in a δ¹⁸O - δD plot (Figure 12), where the Global Meteoric Water Line (GMWL) is also represented. Isotopic data from waters of three nearby rivers Ures (UR), Aconchi (AR) and Sonora (SR) are also plotted as reference of the local isotopic composition of cold water in the study area. As seen in Figure 12, the isotopic compositions of some hot spring waters lie near to the GMWL, particularly those corresponding to the hot springs TCH and MTP from the S and C hydrothermal zones, respectively. The remaining samples (TNB, BCD, ACH, GRN, DVS, HSB, CMP, SMR, TCP, and ARV) show a clear shift for both δ¹⁸O and δD, which is described by the dashed mixing line (Figure 12). This signature is generally the result of mixing with isotopically heavier water. The hot spring waters are probably a mixture of waters similar to the most depleted local groundwaters with a composition that results from the intersection of the mixing line with the GMWL, and an isotopically heavier water not identified in the study zone.

Figure 12 Plot of δD versus δ¹⁸O of all hot spring water samples collected in the four hydrothermal zones (NW, NE, C and S). AR, UR and SR are water samples collected in nearby rivers of Sonora. 

Fluid geothermometry

Pseudo-equilibrium temperature estimates were obtained by the application of seven solute geothermometers (Na/K: 3, K-Mg, Na-K-Ca, Na-Li, and SiO₂). The results of these calculations are reported in Table 7.

Na/K geothermometer

A good agreement among the temperatures calculated by three different versions of the Na/K geothermometer was roughly obtained. The Na/K ratios ranged between 20.7 and 48.1 for the spring waters of the NW, NE (except BCD: 143 and HSB: 467), and S zones, which provided reliable temperature approaches for the reservoir. These estimates must be considered as lower limits of the reservoir temperatures because the chemical signatures of the hot spring waters show evidence of mixing with shallow groundwater, as well as the interaction with surface precipitated minerals. In this context, hot spring waters that are subject to conductive cooling, lateral flow or near surface reactions tend to exhibit high Na/K ratios (between 12.3 and 60) for temperatures ranging from 100 °C to 200 °C (^{Nicholson, 1993}; ^{O’Brien, 2010}). It is well-known that the prediction capability of this geothermometer at temperatures below to 200 °C (i.e., with Na/K ratios ranging from 13 to 27) provide greater uncertainties than those ratios (ranging between 3 and 13) corresponding to high temperatures (between 200 °C and 340 °C). This limitation of the Na/K geothermometer is due to the long WRI times needed for achieving equilibrium between Na and K in controlled laboratory experiments at low-to-medium temperatures (<200 °C; ^{Pérez-Zarate et al., 2015}).

The TNaKVS97 geothermometer systematically provided higher temperature estimates, whereas lower temperature values were computed with the TNaKDSR08 geothermometer. Anomalous higher temperatures of 328±62 °C and 462±87 °C were also calculated for the ARV and MTP hot spring, respectively (values that were excluded from the Table 7). These unrealistic temperature estimates were expected because such samples exhibited anomalous chemical and isotopic signatures (probably caused by the high content of Mg), in comparison with the remaining water samples. On the other hand, temperatures corresponding to the HSB and BCD (which are alkaline together with TCP) hot springs could not be determined with confidence due to the complex chemical and isotopic signatures that revealed some mixing or diluting effects caused by the local hydrogeochemical processes recorded by these samples (see Figures 8-12).

To analyse the anomalous temperature estimates, the Na/K concentration ratios of the geothermal fluids were computed. The highest values of the Na/K ratios were obtained for the HSB and BCD (466.7 and 142.9, respectively), whereas the lowest values correspond to the MTP and ARV hot springs (1.4 and 3.4, respectively). High Na/K ratios may suggest either a preferential leaching of Na or removal of some of the K in the form of secondary alteration products (^{Giggenbach, 1988}), which could be confirmed in future prospection studies by collecting fresh and altered rock samples for analysing the presence of alteration minerals (e.g., muscovite, microcline). Low Na/K ratios were associated to a prediction of unrealistic high temperatures. These geochemical anomalies may be attributed either to ion exchange processes among clay minerals or simply to the lack of equilibrium conditions between solutes and alteration minerals (which is evidenced by the presence of incompletely equilibrated spring waters, probably due to mixing or dilution processes); see the classical Na-K-Mg geothermometer diagram (Figure 13). In this triangular plot, a fast equilibrating process is represented by the K/Mg geothermometer, whereas the slow equilibrating process is given by the Na/K geothermometer (^{Giggenbach 1988}).

Figure 13 Ternary equilibrium diagram showing the relative Na, K and Mg contents for all the hot spring samples collected in the four hydrothermal zones of the Sonora geothermal system (based on ^{Giggenbach, 1988}). 

By excluding the extreme temperatures inferred from ARV and MTP (the highest temperature estimates), and from HSB, and BCD thermal springs (the lowest temperature estimates) (see Table 7), and after applying a geochemometric analysis based on outliers detection/rejection under the assumption of a normal distribution, the mean reservoir temperature predicted with the TNaKVS97 equation for the hot springs CMP, GRN, SMR, TCP, TCH, DVS, TNB, and ACH was 142±16 °C (see Figure 14).

Figure 14 Results of the statistical outlier detection/rejection analysis applied for all the solute geothermometers, and mean reservoir temperatures estimated by each solute geothermometer. 

For the outlier detection/rejection, the univariate analysis applied the Dixon, Grubbs, skewness and kurtosis statistical tests at the 99% confidence level by using the computer code UDASYS of ^{Verma et al. (2013)}. A distribution map showing the most attractive geothermal potential zones (NW and NE) based on these reservoir temperature approaches is shown in Figure 15a. According to this map, lower reservoir temperatures are mainly distributed in the S hydrothermal zone (SMR, TCP and TCH).

Figure 15 Distribution map showing the most attractive geothermal zones (NW and NE) based on the most reliable reservoir temperature approaches (inferred from the solute geothermometers): a) TNaKVS97; b) TNaLiVS97; and c) TNaKCaFT73. The distribution of the surface temperatures measured in the four hydrothermal zones is also plotted as reference (d). 

Mixing geochemical model

Taking into consideration the water mixing processes evidenced from the classification diagrams (Figures 8, 9, 11) and the D/H and ¹⁸O/¹⁶O plot (Figure 12), a mixing geochemical model based on a silica and fluid enthalpy plot was applied for inferring the water composition and temperature at the original reservoir conditions. The silica content of ascending hot spring waters (which fall outside of the full equilibrium curve: Figure 13) and their corresponding fluid enthalpy were plotted in Figure 16. The influence of a conductive cooling was quantitatively examined for determining the mixing degree with colder waters.

Figure 16 SiO₂-enthalpy mixing model and mixing line obtained for the hot spring waters of the Sonora geothermal system. 

After applying the model, it was inferred that the original reservoir conditions indicate a deep reservoir temperature of 210±11 °C, together with a silica content of 300±27 mg/kg, and a fluid enthalpy of 950±85 kJ/kg (see the blue point C in Figure 16). This reservoir temperature is higher than those estimates predicted from the solute geothermometers (Table 7), and it may be considered as the maximum reservoir temperature of the Sonora geothermal system.

According to the Li-Rb-Cs diagram (Figure 11), this reservoir temperature of 210±11 °C is in good agreement with the temperature that controls the illite mineral interaction, and the temperature interval inferred from the K/Rb and K/Cs ratios. All the geochemical signatures observed at deep reservoir conditions reflect the degree of mixing of the waters during their ascent towards the surface.

Fluid-mineral equilibria

For evaluating the pseudo-equilibrium state that actually exhibit most of the waters sampled in the twelve hot springs, the Na-K-Mg diagram was again examined (Figure 13). This diagram associates the hydrothermal fluids with their most probable reservoir temperatures. As seen in Figure 13, most of the hot spring waters fall in the region of mixed waters or partial equilibrium (e.g., ACH and CMP from the NW zone; BCD, GRN, TNB, and DVS from the NE zone; and TCP, SMR and TCH from the S zone). Some exceptions were given by water samples falling outside or in the limits of the full equilibrium line (HSB and BCD from the NE zone), and two remaining water samples located in the section of immature waters (MTP and ARV from the C zone). This geochemical shift was probably caused by the loss of some ions (e.g., K) by these two water samples. Another possible cause may be related to the total dissolved carbonate concentration, which in the specific case of the HSB and BCD water samples comprises both HCO₃ ^- and CO₃ ^2- species). The total carbonate concentration may also be influenced by the CO₂ partial pressure in the deep geothermal fluid (P _CO2) and the pH (^{Nicholson, 1993}). The loss of CO₂ during boiling processes normally raises the pH of the fluid by consuming protons through the involved reactions (i.e., the water usually becomes more alkaline, as in the case of HSB and BCD samples which exhibited the highest pH values of 9.0 and 8.5, respectively; see Table 5).

The lack of equilibrium with the surrounding host rocks is also confirmed by the mixing model (see Figure 16), and explained either as a result of a mixing between hydrothermal waters and colder shallow fluids during the path towards the surface, or simply due to meteoric waters coming from different sources or river waters (Figure 11; e.g., MTP and ARV). The origin of the HSB hot spring water (located above the full equilibrium line in Figure 13) is probably related to mixing with some evaporated water sources (Figure 11).

From the thermal point of view, most of the water samples were distributed in the partial equilibrium region where the reservoir temperatures vary from 75 °C to 175 °C, most of them falling close to the 150 °C isotherm, which seems to be the lower limit of the reservoir temperatures of the geothermal system.

Using 150 °C and 200 °C as rounded-off reservoir temperatures (or min-max values) and the Geochemist’s Workbench software, geochemical equilibrium modelling based on fluid-mineral stability diagrams was carried out. The dominant species and mineral phases that are stable over this temperature interval (i.e., mineral stability diagrams of the primary and secondary minerals for these reactions) were investigated. The following water-rock (mineral) systems were found as the major chemical and mineral signatures underlying the geothermal system under study (see the mineral stability diagrams in Figures 17a and 17b):

Figure 17 Fluid-mineral stability diagrams obtained for the (a) Na₂O-K₂O-Al₂O₃-SiO₂-H₂O and (b) CaO-K₂O-Al₂O₃-SiO₂-H₂O systems. 

1) A Na₂O-K₂O-Al₂O₃-SiO₂-H₂O system, which exhibited ionic exchange between Na and K that seems to be controlled either by the albite/muscovite or the albite/K-feldspar mineral interactions under metastable equilibrium conditions. This system reflects the hydrothermal alteration of primary plagioclase minerals in agreement with the following reaction:

The discharge fluid composition tends to be placed along an extension from either the albite/K-feldspar limit through the muscovite field, or through phases chemically similar like sericite and illite (^{Hedenquist, 1991}). This behaviour can be interpreted as the result of the primary dissolution of feldspars in the reservoir under metastable equilibrium conditions between albite and K-feldspar (^{Torres-Alvarado, 2002}). As seen in Figure 17, the albite and muscovite (K-mica as illite or sericite) interaction shows much better equilibrium conditions at temperatures close to 200 °C than those data corresponding to 150 °C (e.g., ^{Henley and Ellis, 1983}; ^{Hedenquist et al., 2000}). These interactions were mainly observed in the hot spring waters of the NW (ACH and CMP), NE (BCD, GRN and TNB), and S (SMR and TCH) zones.

2) A CaO-K₂O-Al₂O₃-SiO₂-H₂O system, where the ionic exchange between Ca and K is described by the prehnite-laumontite-muscovite mineral assemblage (see Figure 17). These WRI processes seem to be associated with the hot spring waters from the NE (DVS and HSB), and S (TCP) zones. Laumontite is a Ca-rich zeolite typically observed in some geothermal fields at temperatures between 170 °C and 250 °C (^{Reyes, 1990}; ^{Torres-Rodríguez, 2000}). It would be expected that hot spring waters in equilibrium with laumontite indicate the existence of a reservoir with a temperature interval between 120 °C and 210 °C (^{Reyes, 1990}).

Most of these minerals belong to mineral assemblages of low grade metamorphism (^{Liou et al., 1985}). This WRI signature seems to indicate that, according to the regional geology of the study area, the hot spring waters are reaching the chemical equilibria at the contact zones between the plutonic rocks (from the Laramide Batholith) and the tectonic basins (i.e., at regional fault zones under the influence of those plutonic rocks and low-grade metamorphism rocks). The lithological units that represent the former metamorphic features are the intermediate to acidic volcanic rocks from the Tarahumara Formation. The interpretation of the most probable pseudo- or equilibrium conditions among the hot spring waters and the mineral assemblages in the Na₂O-K₂O-Al₂O₃-SiO₂-H₂O and CaO-K₂O-Al₂O₃-SiO₂-H₂O systems seem to be in a good agreement with the expected reservoir temperatures observed in the hydrothermal zones, and the theoretical thermal stability of hydrothermal minerals by ^{Henley and Ellis (1983)} and ^{Hedenquist et al. (2000)}: Figure 18.

Figure 18 Thermal stability of various hydrothermal minerals (modified after ^{Hedenquist et al., 2000} and ^{Henley and Ellis, 1983}). (a) Mineral assemblages associated to pH and temperature; (b) Mineral assemblages associated only to temperature. Dashed orange lines indicate the pH and estimated temperature values prevailing in the hot spring waters reported this study (purple and red dashed lines correspond to minimum and maximum estimates respectively).