This study conducted zircon dating of granite samples HSD-1, HSD-4, HSD-8, HSD-14, HZ-1, and HZ-3. These six samples (∼5 kg each) were processed at the Hebei Institute of Regional Geology and Mineral Resources Survey.

Six zircon grains were separated from these samples using conventional heavy-liquid and magnetic separation techniques. Then, about 100 zircon grains with the highest quality were hand-picked under a binocular microscope. These grains were mounted in epoxy and then polished. Afterward, the grains were photographed using optical microscopy and cathodoluminescence (CL) imaging to reveal their internal morphologies, which were used to select grains and analytical spots. CL images were obtained using a HITACHI S3000-N scanning electron microscope equipped with a Robinson backscattered electron detector and a Gatan Chroma CL imaging system.

To determine zircon U-Pb ages, these above samples were analyzed using an Agilent 7,900 quadrupole ICP-MS equipped with a 193 nm coherent Ar-F laser and a Resonetics S155 ablation cell at the Tianjin Institute of Geology and Mineral Resources. NIST SRM 610 glass used for Pb correction was analyzed after every 15 test samples. The mass bias of Th/U and Pb/Th ratios caused by downhole fractionation and instrumental drift was corrected using zircon reference material 91,500 following instructions of Wiedenbeck et al. (1995). Each zircon analysis involved 15 s blank gas measurement and 50 s analysis. The analytical spots were 29 μm each in size, and laser was emitted at a frequency of 7 Hz, yielding an energy density of approximately 3 J/cm². The grains ablated by laser were carried by He gas into the chamber at a flow rate of 0.35 L/min, where they were mixed with argon gas before being transported to the plasma torch. Zircon standards Temora (Black et al., 2003) and Plesovice (Sláma et al., 2008) were employed in this study. Data processing was performed using Concordia intercept ages on the Tera-Wasserburg plot which was prepared using ISOPLOT (Ludwig, v. 3.75, 2012). Data reduction was conducted using ICPMSDataCal (Liu et al., 2010). The random and systematic uncertainties were estimated using the method proposed by Paton et al. (2010).

Based on LA-ICP-MS zircon U-Pb dating, the cathodoluminescence (CL) diagram of zircon was plotted. The laser denudation was conducted under a beam spot diameter of 40 μm and a duration of 30 s, with GJ-1 used as the external reference material for the calculation of Hf (Geng et al., 2011). In the solution injection mode, a 200×10⁻⁹ standard solution was used. It was fed into a nebulizer at a flow rate of 50 μL/min using a peristaltic pump. Nine Faraday cups simultaneously received ¹⁷²Yb, ¹⁷³Yb, ¹⁷⁵Lu, ¹⁷⁶(Yb+Hf+Lu), ¹⁷⁷Hf, ¹⁷⁸Hf, ¹⁷⁹Hf, ¹⁸⁰Hf, and ¹⁸²W. Data acquisition was performed using the virtual amplifier technology of NEPTUNE (MC-ICPMS). The software automatically replaced the amplifier circuits of the analyzer after the acquisition of a set of data. After nine sets of data were acquired, the amplifier circuits were identical to those of the original analyzer. This technique can be used to effectively eliminate the isotope ratio error caused by the difference in the gain of each Faraday cup receiver, thus improving the accuracy of isotope ratio determination.

The laser ablation test began with the targeting, polishing, and photographing of the staghornite to be tested. Samples were arranged on the double-sided adhesive of the slide and put on the PVC ring. Then, epoxy resin and hardener were mixed thoroughly and injected into the PVC ring. After the resin was fully cured, sample holders were peeled off from the slide and then polished, followed by the taking of the cathodoluminescence, reflected light, and transmitted light photographs of the samples on the target. Based on these three kinds of photographs, the appropriate area of the struvite was selected, and the struvite was stripped using a 193 nm FX laser with a spot diameter of 35, 50, or 75 (um), an energy density of 10–11 J/cm², and a frequency of 8–10 Hz. The laser-stripped materials were fed into the Neptune (MC-ICPMS) with high-purity He as the carrier gas. The receiver configuration was the same as that of the solution injection.

Seven fresh samples were cut into blocks, each measuring ∼70×50×30 mm in size. These blocks were washed with tap water and ground to <200 meshes using an agate mill. Both major and trace element analyses were performed using the method proposed by Liu et al. (2016) at the Central Laboratory of China Railway Resources Group. The contents of major elements were measured using an X-ray fluorescence (XRF) spectrometer (Thermo Fisher ARL Advant'X), with the loss on ignition (LOI) being determined with an electronic analytical balance (CPA225D). The FeO content in the samples was measured using conventional wet chemical titration. The contents of trace elements, including rare earth elements (REEs), were determined using ICP-MS (Thermo Fisher X-series 2) after complete digestion. A detailed account of the test procedures and the lower detection limit were documented in Chinese national standard GB/T14506-2010 (Lv, 2010). The analytical errors varied in a range of 1%–3% of the present values. The whole-rock geochemical analyses were conducted under temperatures of 20–27°C and humidity of 30%–58%.

As presented in Figure 5, zircons from sample HZ-1 were typically transparent, ranging from colorless to slightly brown, and exhibited rectangular to prismatic crystal structures. They were 80–180 μm long, with aspect ratios ranging from 1.5:1 to 3:1. Furthermore, these crystals commonly exhibited oscillatory zoning (Figure 5). Zircons from sample HZ-3 were predominantly transparent, clear to pale yellow, and euhedral to subhedral crystals. They were 70–170 μm long, with aspect ratios varying in the range of 2:1–3:1. Among them, euhedral grains showed concentric zoning with relatively bright cores, as shown in CL images (Figure 5). Zircons from samples HSD-1 and HSD-8 were similar in shape and color. They were slightly larger (typical length: 100–200 μm), with aspect ratios between 1:1 and 3:1, and exhibited weakly oscillatory zoning. Furthermore, zircons from sample HSD-4 were considerably smaller (length: 40–100 μm), with aspect ratios of 1:1–4:1.

The LA-ICP-MS zircon U-Pb ages of the HGF are shown in Figure 6; Table 2. Fifteen analyses showed that sample HZ-1 had a²³⁸U/²⁰⁶Pb age of 135 ± 4 Ma (MSWD = 5.7), while 20 analyses revealed that sample HZ-3 had a weighted average ²³⁸U/²⁰⁶Pb age of 135.3 ± 2.4 Ma (MSWD = 3.2). The two ages, which were consistent within the given error range, indicate that both samples formed during the Cretaceous. In contrast, samples HSD-1, HSD-4, and HSD-8 exhibited significantly older weighted average ²³⁸U/²⁰⁶Pb ages of 152.7 ± 2.7 Ma (MSWD = 4.8), 155.8 ± 3.3 Ma (MSWD = 5.2), and 153 ± 3 Ma (MSWD = 3.8), respectively. These ages, which were also consistent within the given error range, suggest that these three samples intruded during the Jurassic. Sample HSD-14 was emplaced during the Indosinian, as indicated by four data points yielding a²³⁸U/²⁰⁶Pb age of 153.7 ± 1.2 (MSWD = 2.2; probability = 0.002). As shown in Figure 7; Table 3, sample HSD-14 exhibited the highest Ti-in-zircon temperatures from 684°C to 781°C (Ferry and Watson, 2007), with an average of 742°C, while the other samples had lower Ti-in-zircon temperatures, with averages ranging from 617°C to 687°C.

We analyzed the major and trace element compositions of seven representative samples (Table 5). These samples featured high SiO₂ (74.4–76.56 wt.%) and K₂O (3.76–4.84 wt.%) contents. All the samples fell within the high-K calc-alkaline zone (Figure 8A). Samples HSD-4, HSD-8, and HZ-1 exhibited aluminum saturation index (ASI) values of 1.02, 1.05, and 1.05, respectively, all of which were less than 1.1 (Figure 8B). In comparison, samples HSD-1, HZ-3-1, and HZ-3-2 had ASI values of 1.13, 1.13, and 1.11, respectively, all of which were greater than 1.1.

Compared to the other samples, samples HSD-1, HSD-8, and HSD-14 had higher K₂O, CaO, MgO, TiO₂, Zr, Sr, and Ba contents but lower Na₂O and Al₂O₃ contents (Figures 8–11; Table 5). All these samples were depleted in P₂O₅ (Figure 11).

As shown in the primitive mantle-normalized diagrams (Figure 12; Table 5), sample HSD-14 exhibited the strongest depletion in Ti, Sr, La, and Ce but notable enrichment in Rb, U, and Nd. Compared with samples HSD-1 and HSD-8, samples HSD-4, HZ-1, and HZ-3 were relatively depleted in Ti, Sr, and Ba but enriched in Rb, U, Nb, and Nd.

The REE patterns (Figure 13; Table 5) showed that sample HSD-14 exhibited a leftward REE pattern with the lowest La_N/Yb_N ratio (0.56); samples HZ-3-1 and HZ-3-2 showed relatively flat patterns with similar La_N/Yb_N ratios of 1.47 and 1.32, respectively; sample HZ-1 had a much higher La_N/Yb_N ratio of 9.7; among the remaining samples, HSD-1 displayed the highest La_N/Yb_N ratio (5.06), followed by HSD-8 (2.95), while HSD-4 had the lowest La_N/Yb_N ratio (1.19); samples HSD-4, HSD-14, HZ-3-1, and HZ-3-2 exhibited the lowest δEu values, which were 0.1, 0.1, 0.09, and 0.07, respectively, while sample HSD-1 had the highest δEu value, which was 0.32. The Indosinian granite (sample HSD-14) exhibited relatively low total REE contents (ΣREE; 78.38). In contrast, the Cretaceous granites showed the highest ΣREEs, with values of 188.55 ppm, 111.12 ppm, and 149.78 ppm for samples HZ-1, HZ-2, and HZ-3, respectively.

This study conducted the Lu-Hf isotope analyses for zircons that were subjected to U-Pb dating, as shown in Table 6; Figure 14. Specifically, Lu-Hf isotope analysis was conducted on 15 and 20 zircon grains from HZ-1 and HZ-3, respectively. These zircons show a wide range of¹⁷⁶Yb/¹⁷⁷Hf ratios from 0.017479 to 0.126740 (average: 0.050829) and ¹⁷⁶Lu/¹⁷⁷Hf ratios between 0.000509 and 0.003495 (average: 0.00132), indicating low content of radiogenic Hf. Furthermore, 35 analytical spots yielded ε_Hf(t) values of −11.4334 to −3.0396, which corresponded to two-stage Hf model ages (T_DM2) of 1,385–1907 Ma dominated by 1,500–1700 Ma, as depicted in Figure 15.

The zircon ages of the granites in the HGF are summarized in Table 4; Figure 16. These ages indicate three stages of granite emplacement: 1) the Indosinian stage, corresponding to Permian granites (251 ± 9.1 to 253 ± 5 Ma); 2) the Yanshanian stage, including Jurassic granites (152.7 ± 2.7–176.7 ± 1.8 Ma); 3) the Cretaceous stage (135 ± 4–143.6 ± 2.8 Ma) (Xiao et al., 2019). These ages indicate intense Mesozoic (Yanshanian) magmatic activity in the study area. The Cretaceous granites dominate the HGF, exhibiting the largest outcrop area (100 km²; Figure 2), while the Permian granites only occur in the northern part of the geothermal field (Xiao et al., 2019). As indicated by the zircon dating results of samples HZ-1 and HZ-3 taken from the deep part of well HR-1, the surrounding rocks of the geothermal reservoirs in the HGF consist of Cretaceous granites, although Permian granites are the closest to geothermal wells (HR-1 and ZK8).

As shown in the Harker diagram, the samples collected in this study were the product of felsic magmas and had undergone significant fractionation. The Yanshanian granites, represented by the samples collected in this study, show a downward trend of Al, Mg, and Ca contents with an increase in the SiO₂ content (degree of fractionation; Figure 9). The depletion in Sr, Ba, and Ti, as well as prominent negative Eu anomalies (Figure 10), is indicative of the fractionation of plagioclase and Ti-Fe oxides. Figure 11 illustrates decreases in the Zr/Hf, Nb/Ta, and Th/U ratios but increases in the Rb/Sr ratio with an increase in the degree of magmatic fractionation (Stepanov et al., 2016; Breiter and Škoda, 2017; Cai et al., 2020). Furthermore, samples collected from or around geothermal wells exhibited significantly V-shaped REE patterns than those collected at a greater distance from geothermal wells, indicating that they underwent a higher degree of crystal fractionation.

Generally, granites can be categorized into S-, I-, A-, and M-types based on their genesis (Collins et al., 1982; Chappell and White, 2001; Chappell et al., 2012; Zhou et al., 2021). S-type granites originate from sedimentary rocks (Chappell et al., 1999), while I-type granites originate from ancient igneous rocks (Chappell and White, 2001). A-type granites are generally formed by the partial melting of dry granulitic residues in the lower crust, featuring a high abundance of large highly charged cations (Nb, Ga, and REEs) and low Mg, Ca, and Al content (Collins et al., 1982). M-type granites are derived from the subducted oceanic crust or overlying mantle (Whalen et al., 1987). A/CNK serves as an effective factor used to distinguish between S- and I-type granites (Chappell and White, 2001). Sample HSD-14 (Permian granites) had an A/CNK value of 1.2, suggesting that it was S-type granite. This result is consistent with the presence of primary muscovite (Figure 4). In contrast, Yanshanian granites (Jurassic and Cretaceous stages) exhibit A/CNK values of 1.05–1.13, ruling out the possibility of them being S-type granites. As presented in Figure 17, these Yanshanian granites fall within the zone of fractionated granites, suggesting that these granites are highly fractionated I-type granites rather than A-type granites. This finding is consistent with the downward trend of P content with an increase in the SiO₂ content (Figure 10) and can also be evidenced by low Ti-in-zircon temperatures of 475–781°C (Table 3; Figure 7).

Zircon Lu-Hf isotope compositions are an important tool for identifying the magma source of granites (Bhattacharya and Janwari, 2015; Bea et al., 2018). As indicated by the Hf isotope composition and T_DM2 ages (1,385–1907 Ma) of samples HZ-1 and HZ-3, collected at depths of 2,637 m and 2,702 m, respectively (Table 6; Figures 14, 15), the reservoir granites originated from the mixing of Meso-to Paleo-Proterozoic lower crust and juvenile mantle materials. Two theories have been proposed for the formation of highly fractionated I-type granites, namely, the partial melting of crustal materials (Chappell et al., 1999; Wyborn et al., 2001) and the fractional crystallization of mafic melts (Chappell et al., 1999; Chappell et al., 2012). The latter assumption can be excluded for the following reasons: 1) such a large volume of granitoids (80%) is difficult to form due to the too small quantity of mafic rocks within igneous rocks (Zhang et al., 2015); 2) the parent magmas (peraluminous melts) of I-type granites are unlikely to be produced through the fractional removal of hornblende (Chappell et al., 2012). Instead, Cretaceous I-type granites in the study area are likely derived from the partial melting of the lower crust. Since the Early Jurassic (∼190 Ma), the Cathaysia Block had experienced asthenospheric upwelling and intra-continental lithosphere extension after the Indosinian flat-slab subduction of the Paleo-Pacific plate (Li and Li, 2007; Chen et al., 2008). As presented in Figure 18, these granites fall within the plate granite (WPG) zone, indicating an extensional setting. This finding is consistent with the Jurassic-Cretaceous coastward migration of extensional magmatism in the South China Block (Li and Li, 2007; Qiu et al., 2017), during which mafic magmas produced by the partial melting of the asthenospheric mantle triggered the partial melting of lower crust materials and the following generation of felsic magmas. This conclusion accords with the presence of gabbro-basalt rocks (height: 5 km) in the middle crust (depth: 20 km) (Zhang et al., 2005; Zhang et al., 2013).

Granites that have heat production rates >5 μW/m³ are classified as high-heat-producing granites, which have significant effects on the surface heat flow. The presence of high-heat-producing granites at depths of 1,650–2,980 m can increase the surface heat flow from the background value of 68 mW/m² to 93 mW/m², with an increase by 40% (Wang et al., 2023). Compared to magmas, radioactive heat causes a much lower increase in temperature but an extremely longer timespan of thermal anomalies. Mclaren et al. (1999) attributed the temperature increase of 40°C in the contact zone between plutons and host rocks to the decay of radiogenic elements. In the HGF, the Yanshanian granites exhibit heat production rates of 4.51–7.99 μW/m³, averaging 6.35 μW/m³ (Table 5). Therefore, they are high-heat-producing granites. Furthermore, the Cretaceous and Jurassic granites have average heat production rates of 6.63 μW/m³ and 6.06 μW/m³, respectively (Table 5).

According to lithological records and zircon ages (Tables 2, 6; Figures 6, 16), the concealed Cretaceous granites in the HGF occur at depths of 1,560–3,000 m in geothermal well HR-1 (Li et al., 2020b). The large volume of high-heat-producing Cretaceous granites (geothermal reservoir rocks) is an important heat source of geothermal resources in the study area. The surface heat flow of the HGF is 115.5 mW/m², about 1.6 times the background heat flow (73 mW/m²) (Xi et al., 2021). In sum, the heat production of Cretaceous granites plays a major role in the formation of geothermal resources in the HGF.

Previous studies summarized the regional petrogeological characteristics of granites based on granite samples collected from outcrops at a greater distance from HGF (G1, G2, G3, and G5 in Figure 2; Xiao et al., 2019). By comparison, all the samples used in this study were collected from or around geothermal wells (e.g., HR-1) to figure out the petrogenesis of geothermal reservoir rocks buried at depth and the magmatism evolution model of the HGF. As shown by the sampling locations and zircon ages of the samples (Figures 2, 5; Table 4), the ages of the Yanshanian granites exhibit a downward trend from the periphery to the center of the HGF (geothermal wells HR-1 and ZK8). In other words, the geothermal resources are hosted by relatively younger Yanshanian granites. Most especially, for Cretaceous granites, Xiao et al. (2019) obtained the emplacement ages ranging from 138.4 to 143.6 Ma based on the zircon dating of granite samples from G1, G2, and 300 m at the depth of ZK8 (Figure 2; Table 4). This study discovered that the reservoir granites at depth (∼3,000 m) intruded later than the upper and outer Cretaceous granites, with an emplacement age of ∼135 Ma. The reservoir granites exhibit the flattest REE pattern, indicating the strongest fractional crystallization among Cretaceous granites. This result accords with their higher heat generation. Their average heat production (HZ-1 and HZ-3, Table 4) is 6.44 μW/m³, while that of the upper (Id-33, 3.16 μW/m³) and outer (4.51 μW/m³) granites is much lower (representative samples: G1, G2; Figures 10–13). Regarding the Jurassic granites, two zircon ages were obtained by Xiao et al. (2019): 155.8 Ma for G3 and 176.7 Ma for G5. These ages indicate two episodes of Jurassic magmatism in the HGF. The intrusion ages of the Jurassic granites (HSD-1, HSD-4) exposed around the geothermal wells are consistent with those of granites (HSD-8) collected from G3 within the error range. This suggests that the Jurassic granites around the geothermal wells should be formed during the later Jurassic magmatic event (JS 2 in Figure 16) recorded by Xiao et al. (2019). Their average heat production rate (5.99 μW/m³, average value from HSD-1, HSD-4, HSD-8, and Id-04) is twice that (2.62 μW/m³) of granodiorites formed during the former Jurassic magmatic event. It is noteworthy that fine-grained equigranular biotite granites (HSD-4) intruded into HSD-1 (Figure 3) exhibit the strongest depletion of Ba, Nb, Sr, and Ti (Figure 12) and significant Eu anomalies (Figure 13) in the Jurassic granites, suggesting that they have undergone the strongest fractional differentiation (Chen et al., 2014; Liao et al., 2021). This implies that these granites are the latest phase of Jurassic granites.

In sum, the Yanshanian granites around geothermal wells were formed at the late stage of both the Jurassic and Cretaceous. The geochronological data of the granites in the HGF were integrated with the evolution model for plutonism in this study (Figure 19), yielding findings as follows. The Permian two-mica granites (G4 in Figure 2; Figure 19) intruded into the Paleozoic sandstone at the center of the HGF at 253–251 Ma. During the Jurassic, two magmatic events occurred. At 176.7 Ma, granodiorites were emplaced in the southeast of this study area (Jurassic stage 1), while the intrusion of biotite granites (G3 and G6 in Figure 2; Figure 19) occurred at 155.8–153 Ma (Jurassic stage 2). During the Cretaceous, three stages of magmatism occurred: 1) stage 1: biotite granites (G1) were exposed in the northwest of the HGF, with ages ranging from 143.6 to 140.9 Ma; 2) stage 2: biotite granites (G1) were in intrusive contact with granodiorites (G5) in the southeast of the HGF; 3) stage 3: biotite granite intruded under the two-mica granites and Paleozoic sandstones at 135–135.3 Ma as the reservoir rocks of geothermal resources in the HGF.