The Xifeng geothermal field is situated in the central Guizhou Province, southwestern China (Figure 1A). The geological location is in the Guiyang Complex Tectonic Deformation Zone related to the central Guizhou uplift within the Yangtze Craton in South China (Figure 1A; Ling et al., 2015, 2017; Long et al., 2017; Li et al., 2019, Li Y. Y. et al., 2021). The basement is made of Meso-to Neoproterozoic metamorphosed marine sedimentary rocks. The cover consists of Neoproterozoic (Sinian) to Triassic marine sedimentary rocks and Jurassic terrestrial sedimentary rocks, which underwent multiple Phanerozoic tectono-thermal activities (Li et al., 2019; Li J. et al., 2021).

The stratigraphic succession of the Xifeng geothermal field consists of, from bottom to top, the Sinian Qingshuijiang, Nantuo, Doushantuo, and Dengying formations, lower Cambrian Niutitang, Mingxinsi, Jindingshan, Qingxudong, Gaotai, and Shilengshui formations and Loushanguan Group, Permian Liangshan Qixia, Maokou, Wujiaping, and Changxing formations, and Triassic Yelang and Maocaopu formations (Figure 1C; Ling et al., 2015; Long et al., 2017; Li et al., 2019). The Dengying, Qingxudong, Gaotai, and Shilengshui formations and Loushanguan Group are made of carbonates, dominantly dolomite. The Qingshuijiang and Nantuo formations consist of sandstone intercalated with limestone and tillite, respectively. The Doushantuo Formation is made of claystone. The Niutitang Formation is comprised of shale, claystone, siliceous phosphorite, and siltstone, and the Mingxinsi Formation comprises sandstone, limestone, claystone, and shale. The Jindingshan Formation consists of siltstone, sandstone, and limestone. The Permian rocks consist of shale and mudstone (Liangshan Formation), limestone (Qixia and Changxing Formations), dolomite and limestone (Maokou Formation), and claystone, sandstone, and limestone (Wujiaping Formation). The Triassic rocks are made of limestone, shale, and mudstone (Yelang Formation), and carbonate, i.e., limestone and dolomite (Maocaopu Formation) (Ling et al., 2015; Ji, 2015; Long et al., 2018, 2019).

The central Guizhou region experienced a series of tectonic events from Sinian to late Cambrian, such as the Yunan Movement (forming the central Guizhou uplift; Mei et al., 2005; Ling et al., 2015; Long et al., 2017), Duyun and Guangxi Movements (uplifting the central Guizhou; Yu and Wang, 1995; Long et al., 2018), and Dongwu Movement as well as a late Permian Emeishan mantle-plume eruption event (Zhou et al., 2002; He et al., 2005; Jian et al., 2009). These tectonic activities resulted in widespread folding and faulting in the study area, especially the followed Indosinian orogeny, Yanshanian, and Himalayan movements (Huang, 1945; Wang and Mo, 1995; Carter et al., 2001; Chen, 2005; Reid et al., 2007; Li J. et al., 2021), which controlled the distribution of geothermal resources.

Faults are well developed in the study area due to multi-stage tectonic activities, such as the Emeishan mantle-plume eruption, Yanshanian orogeny, and Himalayan crust uplift movement (Figure 1B). These faults are mainly thin-skinned with high hydraulic conductivities and resulted in greater exploitation potential relative to other geothermal fields in China. One of the most important faults is the compresso-shear strike-slip Baimadong fault, which is the main conduit and connects the deep heat source (He et al., 2014; Li and Shen, 2014). The fault is around 50 km long, ENE-trending and dipping to SE at angles of 75° towards 90°, and assigned to a regional large-scale deep fault (Figure 1B; Zhang et al., 2017; Li et al., 2019, Li J. et al., 2021). In addition, secondary Neoid active faults such as the Shaba, Chaoyang, and Shitoutian faults are developed within the tectonic system and act as heat transporting channels (Figure 1B; Song et al., 2014). Among them, the Shaba fault is NE-trending and dipping to SE at angle of 70°. The Chaoyang fault is parallel to the Shaba fault with a length of ∼20 km, NE-trending and dipping to SE at angles of 50° towards 70°. The Shitoutian fault with a NW orientation dips to NE at angles of 75° (Ji, 2015). The fault kinematics and orientations control the distribution of geothermal resources in the study region.

The Xifeng geothermal field covers an area of 102.7 km² (Luo, 2020). Approximately 250 exploration and production geothermal wells have been drilled in this area since the 1950s, and the deepest depth is 2,500 m (He et al., 2014; Yang et al., 2018). Moreover, there are numerous hot springs within this area, including the famous Xifeng hot spring. The measured temperature ranges from 39 to 42°C in the Baimadong area to 53–56°C in the Xifeng area (Luo, 2020). The geothermal gradient of the study area is 2.5–3.0°C/hm with the heat flow rate varying from 50 to 75 mW/m² (Song et al., 2014; Wang et al., 2019a). The Dengying Formation is the main reservoir and the upper Jindingshan, Mingxinsi, and Niutitang formations form the cap of the geothermal system (Figure 2; He et al., 2014; Li and Shen, 2014; Song et al., 2014; Cai et al., 2018).

Fieldwork was carried out in the Xifeng area for collecting samples and obtaining relevant field data in April 2021. A total of 16 water samples (six hot spring samples, two geothermal well samples, six river samples, and two cold spring samples) at Xifeng were collected for major ions, REE, and D-O isotope analyses. Sample locations are shown in Figure 1B.

Time sensitive parameters were tested on site using a portable water quality analyzer. Samples were stored in new 500 ml polyethylene bottles that were rinsed with deionized water twice before sampling. All mentioned hydrochemical analyses of the water samples were performed in the Laboratory of Beijing Research Institute of Uranium Geology. Samples used for analysis of cations were acidified after collection through adding Suprapur HNO₃ to bring the pH to below 2. Analysis of major anions and cations was conducted by using Dionex ICS 1100 ion chromatography through Dionex ionpac AS-19HC and CS12A (4 mm × 250 mm) columns, respectively. The analysis of rare earth elements in water samples was documented using a Thermo Scientific ELEMENT XR inductively coupled plasma mass spectrometer (ICP-MS). The instrument was externally calibrated using a multielement standard solution before ICP-MS analysis. The analytical precision was better than 10% for duplicate analysis of the samples. The composition of deuterium (D) and oxygen (¹⁸O) isotopes of collected water samples was analyzed by a MAT 253 mass spectrometer in a continuous flow mode using a Gas-bench II preparation and introduction system. Isotopic data are expressed in the delta (δ) notation as the per mil (‰) deviation relative to the Standard Mean Ocean Water (SMOW); the analytic precisions (1σ) are ±1.0 and ±0.1‰ for δD and δ¹⁸O, respectively.

Chemical geothermometers are helpful to estimate the reservoir equilibrium temperature of the geothermal system by using the distribution and relative contents of various chemical indicators (i.e., dissolved silica, cation, gas, and isotopes; Das et al., 2021; Saibi et al., 2021). Among them, temperature-dependent solubility (silica) and ion-exchange reactions (i.e., Na-K, Na-K-Ca) are the most applied (Das et al., 2021).

Various geothermometers can obtain different estimated reservoir temperatures as each geothermometer documents the last equilibrium of a specific chemical element and is directly affected by processes of boiling, dilution, and precipitation. Cation geothermometers (Na-K, K-Mg, and Na-K-Ca) and silica geothermometers (quartz no steam loss, quartz maximum steam loss, chalcedony no steam loss, chalcedony maximum steam loss, α-cristobalite, and β-cristobalite) applied in this study are shown in Table 1.

The physico-chemical compositions of samples in this study are shown in Table 2. All samples showed alkaline pH values. In almost all the samples, the anions were dominated by HCO₃ ⁻ and the order of abundance was: HCO₃ ^- > SO₄ ^2- > Cl⁻. Among the cations, the main ion was Ca²⁺, and the following order of abundance was Ca²⁺ > Mg²⁺ > Na⁺ > K⁺. As shown in the Piper (1944) diagram (Figure 3), the waters were mostly the Ca-Mg-HCO₃ type except the Xifeng hot spring water (Ca-Mg-HCO₃-SO₄). Waters had a similar varied pattern for concentrations of the cations and anions, except the high SO₄ contents (88.3–91.5 mg/L, avg. = 89.9 mg/L) in the Xifeng hot spring as illustrated in the Schoeller (1995) semi-logarithmic diagram (Figure 4). The ternary plot of Cl-SO₄-HCO₃ (Giggenbach, 1991) was conducted to evaluate the type of water mixed in thermal or non-thermal fluids, notably the peripheral waters involved are shown in Figure 5. The groundwater reacted with host rocks though the unequilibrated state (Figure 6). The intercorrelation of ions during the process of the water-rock reaction is shown in Figure 7.

The REE signatures of geothermal fluids can be used to assess the influences of water-rock interaction. REE composition of sampled waters were analyzed in this study, and the results are shown in Table 2. Four types of water appeared to have similar PAAS-normalized REE patterns, which were featured by heavy rare earth elements (HREE) enrichment compared to light rare earth elements (LREE), and positive Eu anomalies (Figure 8). The total rare earth elements (REE) contents in the geothermal well fluids (∼0.048 to ∼ 0.057 μg/L), hot spring waters (∼0.039–∼0.214 μg/L), river waters (∼0.067 to ∼0.214 μg/L), and cold spring waters (∼0.157 to 0.176 μg/L) were relatively low. The δCe values were 0.346–0.377, 0.133–0.577, 0.619–0.238, and 0.864–1.167, and δEu values were 16.312–28.836, 3.845–44.950, 3.918–19.978, and 2.106–3.723 for waters of the geothermal well, hot spring, river, and cold spring, respectively.

Hydrogen and oxygen isotopic compositions of sampled waters in this study are shown in Table 2. The measured δD_v-SMOW values ranged from −65.6 to −65.6‰ (avg. = −65.6‰, n = 2), −59.1 to −56.2‰ (avg. = −57.6‰, n = 6), −55.6 to −47.9‰ (avg. = −51.8‰, n = 6), and −47.7 to −45.9‰ (avg. = −46.8‰, n = 2) for the geothermal well, hot springs, river, and cold spring waters, respectively. The measured δ ¹⁸O_v-SMOW values were −9.2 to −8.7 (avg. = −8.95‰, n = 2), −8.5 to −7.8‰ (avg. = −8.13‰, n = 6), −8.2 to −7‰ (avg. = −7.71‰, n = 6), and −7 to −6.5‰ (avg. = −6.75‰, n = 2) for the geothermal well, hot spring, river, and cold spring waters, respectively.

Since most sampled waters were assigned to the Ca-Mg-HCO₃ type (Figure 3), this bicarbonate and Ca-Mg dominated feature can be attributed to the interaction with reservoir rocks, which are mainly Sinian dolomites. Ca²⁺, Mg²⁺, and HCO₃ ⁻ were sourced from the dissolution of dolomites following the equation of CaMg ( C O 3 ) 2 ( dolomite ) + 2 H 2 C O 3 = C a 2 + + M g 2 + + 4 HC O 3 − (1)

As illustrated in Figure 6, the water-rock interaction is a dominant process in thermal and cold springs of the Xifeng geothermal field, although water is not equalized with the host reservoir rocks. Water-rock reactions can be further evidenced by the δ¹⁸O results of geothermal waters, which deviate 2‰ of units from LMWL (Figure 11; Taylor, 1977). Additionally, most REE concentrations of sampled waters were above the detection limits (Table 2), indicating that REE in the samples were not derived from meteoric waters (Lewis et al., 1994). The positive Eu anomaly (triggered by the interaction between water and feldspar-bearing rocks or the physicochemical conditions; Sverjensky, 1984; Şener et al., 2017) supports the opinion that REE in the samples were inherited from feldspar-bearing dolomites through the water-rock interaction since the temperatures of geothermal fluids were much lower than 200°C (Figure 8; Sverjensky, 1984).

The characteristics of major ions and their intercorrelation can be used to deduce the geochemical processes caused by the water-rock reaction that the groundwater encounters along its flow path (Adams et al., 2001). According to Eq. 1, the dissolution of dolomite would produce a (Ca²⁺ + Mg²⁺)/HCO₃ ^- molar ratio of 1:2 (Han et al., 2013; Belkhiri and Narany, 2015). However, the ratios of (Ca²⁺ + Mg²⁺)/HCO₃ ⁻ were lower than 0.5 (Figure 7A), indicating that other sources of HCO₃ ^- were involved in the geothermal fluids (e.g., silicate weathering), whereas the deficit in contents of Ca²⁺ and Mg²⁺ should be influenced by other hydrochemical processes (e.g., ion-exchange reaction), rather than the sole dissolution of dolomite. The affection of silicate weathering can also be inferred from the high values of Ca²⁺/Mg²⁺ (>1) in this study (Katz et al., 1997). Moreover, the Na⁺/Cl⁻ molar ratios of the thermal waters were much higher than 1, further confirming the fact that the excess Na⁺ was sourced from silicate weathering (Figure 7B; Bob et al., 2015). The plots for most samples were well above the 1:1 line of (HCO₃ ⁻ + SO₄ ²⁻)/(Ca²⁺ + Mg²⁺) (Figure 7C), implying the occurrence of ion exchange. It is clear that the contents of Ca²⁺ and SO₄ ²⁻ in natural water commonly depended on the dissolution of gypsum and precipitated processes, as shown in the equation of CaSO₄ + 2H₂O = Ca²⁺ + SO₄ ²⁻ + H₂O (Das et al., 2021). However, with influence of silicate weathering, the ratio of Ca²⁺ versus SO₄ ²⁻ was high than 1, except samples collected from the Xifeng hot spring (Ca²⁺/SO₄ ²⁻ < 1; Figure 7D). The high SO₄ ²⁻ in the Xifeng hot spring may be due to the dissolution of other sulphate minerals during the deeper and long-duration circulation from high-elevation SW towards low-elevation NE in the study area (Figure 4; Song et al., 2014).

The Giggenbach Na-K-Mg ternary diagram is used for categorizing waters as full equilibrium, partial equilibrium, or immaturity compared with reservoir rocks (Giggenbach, 1991; Giggenbach and Goguel, 1988). As shown in Figure 9, all the sampled waters were plotted in the field of immature waters, indicating none of them were in full or partial equilibrium with the host rocks or may have mixed with shallow groundwater. It is established that only the waters with features that were fully or partially equilibrated with the host rocks can be used to obtain reliable reservoir temperatures; immature waters produce less reliable results which are usually not regarded as the real reservoir temperatures. Hence, cation geothermometers were not appropriate for the sampled waters in this study. This can be further evidenced by geo-thermometry results obtained through different Na-K geothermometers, which showed higher estimated reservoir temperatures than that of the Na-K-Ca and K-Mg geothermometers (Table 3). The inconsistent results of geothermometers containing Ca and/or Mg with that of Na-K are attributed to the mix of shallow groundwaters and/or surface waters (Figure 5; García-Soto et al., 2016). This mix made the reliability of silica-based geo-thermometry stronger than cation geothermometers (Pandarinath, 2011; Pandarinath and Domínguez-Domínguez, 2015; Saibi et al., 2021). The reliability of silica-based geo-thermometry can be further supported by Figure 10, which shows that the thermal and cold springs of the Xifeng geothermal field were near to the equilibrated curve of chalcedony + quartz (Giggenbach, 1991). Given that the temperatures calculated by chalcedony and cristobalite geothermometers were much lower than the real temperatures of the hot spring and wellhead, reservoir temperatures generated though the quartz (conductive and adiabatic) geothermometers were the most reliable (García-López et al., 2014). The quartz (conductive and adiabatic) geothermometers gave the average estimated reservoir temperature of 77°C. The depth of the geothermal reservoir is around 2.57 km in the Xifeng geothermal field, assuming a temperature gradient of 3°C/100 m (Song et al., 2014).

Oxygen-18 (¹⁸O) and deuterium (D) contents in sampled waters from the Xifeng geothermal field were analyzed to confirm the source and circulation mechanism of geothermal fluids. On the δ¹⁸O-δD diagram (Figure 11), most samples were plotted close to the local meteoric water line (LMWL: δD = 8.83 δ¹⁸O+ 22.15, Zhang et al., 2005) and global meteoric water line (GMWL: δD = 8 δ¹⁸O+ 10, Craig, 1961), indicating a common meteoric water origin. This is in accordance with results of the PAAS-normalized REE patterns (Figure 8) and variable characteristics of ions (Figure 4), which showed a similar source of geothermal fluids.

Craig (1961) established the δD and δ¹⁸O values of precipitation relating to altitude effect, and proposed that meteoric water derived from higher elevation is gradually lighter. In Figure 11, the slightly depleted δD and δ¹⁸O values of geothermal waters compared to those of local river waters testify that the recharge elevation of the former is higher than that of the latter. Generally, the recharge elevation of groundwater can be calculated in terms of the following formula: H = δ 18 O gw − δ 18 O 1 w gra d 18 O + h (2)

In Eq. 2, H (m) = recharge elevation; δ¹⁸O_gw = oxygen isotope value of groundwater; δ¹⁸O_lw = oxygen isotope value of local meteoric water; grad¹⁸O (‰/km) = isotope elevation gradient of meteoric water, and h (m) = elevation of the local meteoric sampling point. The most depleted δ¹⁸O value (−9.2‰) in the geothermal samples was regarded as the δ¹⁸O_gw to minimize the positive isotopic shift effect in this study. The average oxygen isotope value (−7.7‰) of surface meteoric water obtained in this study was used as δ¹⁸O_lw. Hence, the recharge elevation was calculated to be 1,583 m as the δ¹⁸O vertical gradient in Guizhou was assumed to be −3.1‰/km (Yu et al., 1984), and with the elevation of local meteoric sampling point at 1,100 m. Combined with the geological setting, this result is reasonable as the sampling elevation is around 1,250 m for the samples collected from southwest and 750 m for those from northeast and southeast of the geothermal field.

The formation of the Xifeng geothermal field resulted from recharge, deep circulation, and secondary rising meteoric water along the faults. Based on the topographic features and geological conditions (Figures 1, 2), we propose the basic conceptual model for the genesis of the Xifeng geothermal system, as illustrated in Figure 12 and discussed.

Results of stable isotopes and geology revealed that the local meteoric water from the surrounding mountains seemed to infiltrate to depth through faults and crack zones and were the dominant recharge source for the geothermal system. Sedimentary rocks of Lower Cambrian Jindingshan, Mingxinsi, and Niutitang formations overlie the Upper Sinian Dengying reservoir units and act as cap rock (Figure 1C). Moreover, the black shale of the lower Cambrian Niutitang Formation in South China is well known to be enriched in radioactive heat-producing uranium elements, which are 10 and 6–20 times the content compared with that of crustal sedimentary rocks and the crust, respectively (Coveney and Nansheng, 1991; Ni et al., 2012; Pi et al., 2013; Yang et al., 2013; Pagès et al., 2018; Li et al., 2019). At Xifeng, the average U, Th, and K contents for the black shale of the Niutitang Formation is 32.61 (μg/g), 13.25 (μg/g), and 1.86 (%), respectively (Chen 2005; Pi et al., 2013). Generally, the radioactive heat production can be calculated via the formula of A = 0.01ρ (9.52C_U + 2.56C_Th + 3.48C_K) (Rybach, 1976). Where A (μW/m³) refers to radioactive heat production, ρ (g/cm³) is the rock density, C_U (μg/g), C_Th (μg/g), and C_K (%) are the U, Th, and K concentrations, respectively. Thereby, the heat production potential of the lower Cambrian Niutitang Formation in the study area is calculated to be 9.57 μW/m³ with the rock density of 2.728 (Zhao et al., 1995). This value is higher than the average heat production of the continental crust (∼1 μW/m³; Waples, 2001) and can be regarded as an effective radiogenic heat source (Paternoster et al., 2017). Hence, except for the role of cap rock, the Niutitang Formation also served as a heat source due to the producing of radioactive heat (Figure 12). Furthermore, since the multi-period (especially in Yanshanian and Himalayan) active strike-slip of the Baimadong fault and Neoid active faults (e.g., Shaba, Chaoyang, and Shitoutian faults; Figure 1B) are well developed within the tectonic system of the Xifeng geothermal field, the tectonic frictional heat is speculated to act as a supplementary heat source due to the mechanical friction produced along the faults during deformation (Mase and Smith, 1987; He et al., 2014; Song et al., 2014; Zhu, 2016). Moreover, these deep large scale faults also connect to the deep earth and act as a heat transporting channel for deep heat (Figure 12; He et al., 2014; Li and Shen, 2014; Song et al., 2014; Li et al., 2019, Li J. et al., 2021). The existence of acid magmatic intrusion material under the central Guizhou uplift can be considered as the deep heat source, evidenced by the available geophysical information (Zhang et al., 2017).

Once the meteoric water was warmed up, it would interact with the host rocks. The interaction between water and Sinian dolomites is the dominant process for the formation of Ca-Mg-HCO₃ type waters in the Xifeng geothermal field. Moreover, the silicate weathering and ion exchange are also responsible for the formation of geothermal fluids in the Xifeng area. As to the Ca-Mg-HCO₃-SO₄ type in the Xifeng hot spring, it may be due to the dissolution of other sulphate minerals during the deeper, and long-duration circulation from high-elevation SW (i.e., 1,250 m, Baimadong area) towards low-elevation NE (i.e., 750 m, Xifeng area) in the study area (Figure 12; Song et al., 2014). This is consistent with the results of recharge elevation and geological setting.

Then finally, as the deep-infiltrated waters experienced considerable deep circulation, the geothermal fluids rose again along the main ENE-SE, NE-SE, and NW-NE faults. Moreover, the ascending geothermal fluids were mixed with cold groundwater in the subsidiary fractures near the surface (Figure 5; Giggenbach, 1991). The tectonic and stratigraphical features at the study area are favorable for the formation of abundant geothermal resources.