Quantifying Exhumation at the Giant Pulang Porphyry Cu-Au Deposit Using U-Pb-He Dating

The Triassic Pulang porphyry Cu-Au deposit, located in the South Yidun terrane, is the oldest and one of the largest porphyry deposits in the southeastern Tibetan Plateau. The mineralization occurs mostly in the potassic alteration zone of the Pulang intrusive complex. U-Pb-He triple dating, namely apatite (U-Th)/He, zircon U-Pb, and zircon (U-Th)/He dating, together with inverse thermal modeling, reveals that the Pulang complex was emplaced at a paleodepth of ~5.0 to 6.5 km at 215 ± 2 Ma. The deep-level emplacement of the complex, coupled with the episodic replenishment of the magma chamber, gave rise to the establishment of a prolonged magmatic-hydrothermal system at Pulang. Although a range of single-grain zircon and apatite (U-Th)/He ages were obtained on each sample, the weighted mean zircon and apatite (U-Th)/He ages vary systematically with elevation, defining a multistage cooling/denudation history at Pulang. Specifically, three phases of cooling were recognized from inverse thermal modeling, including rapid cooling (80°–120°C/m.y.) in the Late Triassic, moderate cooling (3°–5°C/m.y.) from the Late Triassic to Early Cretaceous, and a protracted slow cooling period (<1°C/m.y.) from the Early Cretaceous to the present day. The first phase of cooling can be mainly attributed to magmatic cooling, whereas the later two phases of cooling were predominantly controlled by uplift and denudation processes. Moreover, the remarkable decrease in the cooling rate from the second to the third phase can be linked to a decreasing erosion rate during the third phase, supported by age-elevation relationships. Overall, our results indicate that the Pulang complex experienced two stages of exhumation at 33 to 45 m/m.y. and 5 to 17 m/m.y. Based on these data, we estimate that approximately 558- to 1,099-m thickness of materials have been removed from the Pulang complex during uplift and erosion, including a large volume of ore. The long time span (>50 m.y.) of extremely slow cooling and erosion at Pulang could be related to the formation and preservation of a peneplain on the southeastern Tibetan Plateau since the Late Cretaceous. A relict peneplain thus signifies a favorable tectonic environment for the preservation of ancient porphyry systems worldwide.


Introduction
Porphyry ore deposits are products of complex magmatic and hydrothermal activity at convergent plate margins (Richards, 2011;Cooke et al., 2014) and commonly form at paleodepths of 1 to 6 km (mode average of ~2 km; Wilkinson and Kesler, 2007;Yanites and Kesler, 2015). Previous compilations of geochronological data revealed that porphyry deposits worldwide typically have Phanerozoic ages, with a particularly strong Cenozoic age population (e.g., Singer et al., 2008). This temporal distribution of porphyry deposits largely reflects the increased likelihood of older deposits having been destroyed by continuous erosion. Quantifying postmineralization exhumation processes has thus become fundamental to understanding how and when hypogene ores were unroofed after formation and when they were exposed at the paleosurface. In addition, the exhumation history of porphyries can provide insight into regional, and even global, tectonic evolution (e.g., Yanites and Kesler, 2015;Zhao et al., 2016). Our understanding of the timing of emplacement and exhumation processes has been facilitated by the emergence of low-temperature thermochronometry methods (e.g., Wolf et al., 1996;Farley et al., 1998).
Fission-track and (U-Th)/He dating are two common thermochronological methods. They have closure temperatures (Tc) ranging from ~60° to 250°C and are thus sensitive to exhumation through crustal depths of one to several kilometers (e.g., Reiners and Brandon, 2006). The apatite (U-Th)/He system has the lowest closure temperature (45°-70°C; Farley et al., 1998;Wolf et al., 1998) and has the potential to place temporal constraints on exhumation processes from 1-to 3-km depth to the surface under a geothermal gradient of 20° to 40°C/km, commensurate with the average depth of most porphyry systems. The apatite (U-Th)/He thermochronometer has been used widely in combination with other thermochronometers to quantify the cooling and exhumation processes at porphyry and/or skarn deposits (e.g., McInnes et al., 1999McInnes et al., , 2005aMasterman et al., 2005;Xu et al., 2006;Harris et al., 2008;Bineli Betsi et al., 2012;Braxton et al., 2012;Zhao et al., 2015Zhao et al., , 2016, orogenic-type Au deposits (Zeng et al., 2013;Liu et al., 2017), and diamondiferous kimberlite occurrences (e.g., McInnes et al., 2009;Evans et al., 2013). McInnes et al. (2005a, b) performed apatite (U-Th)/He, zircon (U-Th)/He, and zircon U-Pb dating on samples from porphyry deposits in Iran, Chile, and Indonesia. Combined with inverse numerical modeling (i.e., 4DTherm; Fu et al., 2010), the data were used to decode thermal histories (magmatic through hydrothermal and exhumation) over a 700°C temperature interval. However, modeling requires knowledge and assumptions of physicochemical parameters of the igneous bodies and their country rocks (e.g., geothermal gradient, size of intrusion, temperature, etc.), which are not always available. The vertical transect method, devised by Wagner and Reimer (1972), has been utilized to deduce simple age-elevation relationships for thermochronometers. This relationship provides a direct basis to estimate long-term erosion rate, independent of geothermal gradient and other parameters, and has been used extensively in the study of geomorphic evolution (e.g., Reiners et al., 2003;Flowers and Farley, 2012).
We have applied the vertical transect sampling strategy, U-Pb-He dating, and numerical modeling to study the cooling and erosion history of the giant Pulang porphyry Cu-Au deposit, which formed during the Late Triassic in the Yidun terrane, at the southeastern margin of the Tibetan Plateau. Pulang was selected for the case study, because it is the oldest and one of the largest porphyry deposits in southeastern Tibet, and several >1,000-m deep drill holes provide excellent vertical transects for sampling. In addition to quantifying the thermal history of the deposit and determining the uplift and erosion rate at Pulang, the research aims to shed light on the tectonic evolution of the southeastern Tibetan Plateau and the preservation of ancient porphyry deposits.

Geology of the Yidun terrane
The Yidun terrane, also known as the Chuanxi Plateau or Litang Plateau, is located between the Qiangtang and Songpan-Ganzi terranes (Fig. 1a). It is characterized by abundant low-relief surfaces (i.e., relict peneplains) at high elevations Tian et al., 2014). The basement of the Yidun terrane is composed of Neoproterozoic granitic gneisses and metavolcanic rocks, which are similar to the lithological assemblages of the Kangding complex in the western margin of the Yangtze block (Fig. 1b), implying that the Yidun terrane was part of the Yangtze block prior to the Early Permian (Chang, 2000;Song et al., 2004). Triassic sequences of clastic rocks, intercalated with arc-type volcanic rocks, overlie Paleozoic sequences (Fig. 1b). These volcanic rocks formed during the westward subduction of the Ganzi-Litang Ocean in the Late Triassic (e.g., Leng et al., 2014). A number of large coeval granite batholiths, and felsic-intermediate porphyry stocks and dikes with associated ore deposits, were emplaced into the upper crust of the Yidun terrane (Reid et al., 2007;Leng et al., 2014;Tian et al., 2014). The absence of Jurassic-Cretaceous strata indicates that this terrane was uplifted above sea level no later than the Early Jurassic.
The Yidun terrane experienced at least two main phases of deformation and metamorphism from Early Triassic to Early Jurassic (Reid et al., 2005;Roger et al., 2011). The early phase of deformation and metamorphism was recorded mainly by Paleozoic sequences in the western Yidun terrane, producing up to amphibolite facies metamorphic rocks. This phase of deformation and metamorphism resulted from northeastward accretion of the eastern Qiangtang terrane onto the Yidun terrane along the southern Jinsha suture zone in the Early Triassic (Reid et al., 2005). Late-phase metamorphic grades were up to greenschist facies and affected Middle-Late Triassic rocks of the eastern Yidun terrane during the Late Triassic to Early Jurassic (Reid et al., 2005;Yang et al., 2012). This stage of metamorphism and deformation was related to the closure of the Ganzi-Litang Ocean and the subsequent collision between the Songpan-Ganzi and Yidun terranes (Reid et al., 2005;Yang et al., 2012).
The Yidun terrane was intruded by the Mesozoic-Cenozoic granitoids that were derived from melting of old continental crust in an intraplate extensional setting (Reid et al., 2007;Wang et al., 2014). However, the emplacement of these Cretaceous felsic plutons did not affect the cooling history of the Triassic plutons (e.g., Reid et al., 2005;Roger et al., 2011). Following collision of the Indian and the Eurasian plates in the Tertiary, the Yidun terrane was incorporated into the modern Tibetan Plateau and deformed by numerous strikeslip fault movements (Wang and Burchfiel, 2000;Zhang et al., 2015Zhang et al., , 2017.

Geology and previous geochronology of the Pulang deposit
The Pulang porphyry deposit, located in the south Yidun terrane, contains proven reserves of approximately 4.31 million tonnes (Mt) of Cu and 113 t of Au, with an average grade of 0.34% and 0.09 g/t, respectively . Detailed geologic descriptions have been provided by Fan and Li (2006) and Li et al. (2011), with a summary presented below.
The Cu-Au orebody is mostly confined within the Pulang intrusive complex, which was emplaced into the Triassic clastic and volcanic rocks of the Tumugou Formation (Fig. 2a). The complex is composed of five phases of porphyry stocks and dikes and covers an area of 16 km 2 . The stocks include premineralization quartz diorite porphyry, synmineralization monzodiorite porphyry and quartz monzonite porphyry stocks, and postmineralization granodiorite porphyry and andesite porphyry dikes. Associated hydrothermal alteration assemblages include potassic, sericite, and propylitic alterations, which were extensively developed in the Pulang complex. Potassic alteration generally occupies a central zone, which is surrounded by propylitic alteration on the periphery. Both alteration types were locally overprinted by late sericite alteration. Ores occur mainly within the potassic alteration zone (Fig. 2a, b). Biotite-magnetite alteration affected adjacent clastic rocks of the Tumugou Formation and produced hornfels zones surrounding the complex (Fig. 2a).
During the past decade, many attempts have been made to date the Pulang complex and ores (Table 1). Wang (2008) obtained zircon SHRIMP U-Pb ages of 228 to 226 Ma for the quartz monzonite porphyry, slightly older than the zircon isotope dilution-thermal ionization mass spectrometry (ID-TIMS) U-Pb age of 211.8 ± 0.5 Ma for the same porphyry (Pang et al., 2009). Pang et al. (2009) also obtained two zircon U-Pb ages of 221.0 ± 1.0 and 206.3 ± 0.7 Ma for the quartz diorite porphyry and granodiorite porphyry, respectively. Similar but less precise U-Pb ages for these porphyries were also provided by Wang et al. (2011) and Liu et al. (2013). Zeng et al. (2006) dated the mineralization, providing a reliable Re-Os age of 213 ± 3.8 Ma for molybdenite. The hornblende and biotite grains separated from some potassic altered porphyries yielded 40 Ar/ 39 Ar plateau ages varying from 216 to 210 Ma (Zeng et al., 2006;Wang, 2008;Li et al., 2011), broadly comparable to the molybdenite Re-Os ages and thus probably recording the timing of ore-related potassic alteration. In contrast, two K-feldspar samples separated from the quartz monzonite porphyry yielded much younger K-Ar and 40 Ar/ 39 Ar plateau ages of ~182 Ma (Zeng et al., 2006;Wang, 2008;Li et al., 2011). Considering that the closure temperature of the K-feldspar Ar-Ar system is as low as ~150° to 350°C (e.g., Foland, 1994;Villa and Hanchar, 2013), these ages likely represent a late thermal event.

Sampling and Analytical Methods
Given that an age vs. elevation profile can be used to estimate exhumation rate independently, we collected seven porphyry samples at approximately 500-m intervals from some    Note: All the 40 Ar/ 39 Ar dates are reported as plateau ages for the K-bearing minerals in this paper, unless otherwise specified Abbreviations: GP = granodiorite porphyry, ID-TIMS = isotope dilution-thermal ionization mass spectrometry, LA-ICP-MS = laser ablation-inductively coupled plasma-mass spectrometry, Moly = molybdenite, QDP = quartz diorite porphyry, QMP = quartz monzonite porphyry, Qtz = quartz representative deep drill holes (e.g., ZK0409 and ZK2812) and one sample from a surface outcrop (Table 2; Fig. 2). Zircon and apatite crystals were separated from crushed rock using a combination of magnetic and heavy liquid separation techniques. Individual crystals were handpicked under a binocular microscope. Zircon and apatite crystals were used for zircon U-Pb and zircon and apatite (U-Th)/He dating. Detailed dating results are provided in Electronic Appendix Table A1. Table 2 gives an overview of sampling locations, lithology, and multiple radiometric ages for the samples.

Zircon LA-ICP-MS U-Pb dating
Zircon U-Pb dating was performed on an Agilent 7500cs quadrupole inductively coupled plasma-mass spectrometer (ICP-MS) with a 193-nm Coherent Ar-F gas laser at the University of Tasmania, Hobart, Australia. Downhole fractionation, instrument drift, and mass bias correction factors for Pb/U ratios were calculated using standard blocks, including two analyses of the primary standard 91500 (Wiedenbeck et al., 1995) and one analysis on each of the secondary standard zircons (Temora, Black et al., 2003;and GJ-1, Jackson et al., 2004), which were analyzed at the beginning of the session and every 15 unknown zircons, using the same spot size and analytical conditions as those used on the samples. Additional secondary standards, Plešovice and Mud Tank, were also analyzed to monitor the analytical accuracy. The correction factor for the 207 Pb/ 206 Pb ratio was calculated using large spots of NIST610 analyzed every 30 unknowns and corrected using the values recommended by Baker et al. (2004). Each analysis began with a 30-s blank gas measurement followed by a further 30 s of analysis with the laser firing. Zircons were sampled using 32-µm spots using a frequency of 5 Hz and an energy density (measured at the sample surface) of approximately 2 J/cm 2 . A flow of He carrier gas at a rate of 0.35 L/min carried ablated material out of the chamber to be mixed with Ar gas and carried to the ICP-MS plasma torch. Isotopes measured were 49 Ti, 56 Fe, 90 Zr, 178 Hf, 202 Hg, 204 Pb, 206 Pb, 207 Pb, 208 Pb, 232 Th, and 238 U, with each element measured every 0.16 s with longer counting times on the Pb isotopes compared to the other elements. Data were reduced using methods described in Meffre et al. (2008) and Sack et al. (2011). Element abundances on zircons were calculated using Zr as the internal standard element, assuming stoichiometric proportions (43.14% Zr in zircon) and using NIST610 to correct for mass bias and drift. Secondary standards reproduced accepted values within 2%. Analytical results for the standards are presented in Electronic Appendix Table A2.

Zircon and apatite (U-Th)/He dating
Zircon and apatite (U-Th)/He dating was conducted at the John de Laeter Centre, Curtin University, Western Australia, using methods detailed in Evans et al. (2005Evans et al. ( , 2013. Zircon and apatite grains were selected by handpicking to avoid grains with cracks or U-and Th-rich mineral inclusions that may contribute excess He to the analysis. Grain measurements were taken for the calculation of an alpha correction factor (Ft; Farley, 2000). Helium was thermally extracted from single crystals of zircon or apatite, loaded into niobium or platinum microcrucibles (respectively), and heated using a diode laser. 4 He abundances were determined by isotope dilution using a pure 3 He spike, calibrated daily against an independent 4 He standard. The uncertainty in the 4 He measurement of sample was <1%. For zircon, the U and Th concentrations were determined using isotope dilution-inductively coupled plasma-mass spectrometry. Zircon grains were removed from the laser chamber and transferred to Parr bombs, where they were spiked with 235 U and 230 Th and digested at 240°C for 40 h in HF. Standard solutions containing the same spike amounts as samples were treated identically, as were a series of unspiked reagent blanks. A second bombing in HCl for 24 h at 200°C ensured dissolution of fluoride salts, and final solutions were diluted to 10% acidity for analysis on an Agilent 7500cs mass spectrometer. For degassed apatite, the U and Th contents were determined by isotope dilution using 235 U and 230 Th spikes. Twenty-five μL of a 50% HNO3 solution containing ~15 ppb of 235 U and 5 ppb of 230 Th was added to each sample. The apatite was digested in the spiked acid for at least 12 h to allow the spike and sample isotopes to equilibrate. Standard solutions containing the same spike amounts as the samples, in addition to 25 μL of a standard solution containing 27.6 ppb U and 28.4 ppb Th, were treated identically, as were a series of unspiked reagent blanks. Prior to analysis on the mass spectrometer, 250 μL of MilliQ water was added.
Both U and Th isotope ratios were measured to a precision of <2%. The (U-Th)/He thermochronology method at Curtin University has an overall precision of <6% for zircon (based on repeated analysis of an in-house standard) and 2.5% for apatite based on multiple age determinations (n = 70) of Durango standard, which produce an average age of 31.5 ±

Inverse thermal modeling
Numerical modeling techniques provide effective and complementary tools for interpreting real thermochronometry data. In this study, two computational modeling software programs (HeFTy and 4DTherm) were independently used to determine the thermal history of the Pulang complex, utilizing dating results determined in this work and previously published data. Here we provide a brief overview of the two approaches with detailed descriptions given in Ketcham (2005) and Fu et al. (2010), respectively. HeFTy is one of the most popular thermal modeling programs for low-temperature thermochronological data, including apatite and zircon (U-Th)/He ages (Ehlers et al., 2005). In the latest version of this program (i.e., HeFTy v.1.9.3), the kinetic models of Flowers et al. (2009) and Guenthner et al. (2013) were used for apatite and zircon (U-Th)/He calculations, respectively. Input data for the inversions included apatite and zircon (U-Th)/He ages, grain dimensions, and U-Th concentrations. Models were constrained by the present-day surface temperature (i.e., 10 ± 2°C) and the time-temperature constraints placed around the published 40 Ar/ 39 Ar ages of hornblende and biotite (Table 1), and the zircon U-Pb ages of this study (Table 2). Models were run until 50 good fits, as defined in Ketcham (2005), were obtained.
4DTherm v.1.2 uses a 2-D explicit finite difference solution that addresses conduction cooling, magma convection, groundwater circulation, latent heat of crystallization, and fusion as well as exhumation and erosion processes (Ehlers et al., 2005;Fu et al., 2010). The main application of this program is to reconstruct cooling and exhumation histories of igneous intrusions, which is achieved by determining the emplacement depth of igneous intrusions and calculating the erosion rates for igneous rocks and corresponding country rocks (Ehlers et al., 2005). The calculation of erosion rates is mainly based on the input apatite and zircon (U-Th)/He data, whereas the emplacement depth and eroded thickness of the igneous body are determined through inverse modeling using an iterative trial-and-error strategy (see Fu et al., 2010). Data inputs to the program include the following: (1) size and shape of igneous and country rock units, (2) residual sizes, position, and properties of igneous bodies, (3) age data and closure temperatures, and (4) sample position. The outputs available after each successful run include digitized agetemperature cooling curves with highlighted key points and a number of physical parameters related to the dynamic processes of magmatic-hydrothermal cooling, such as depth and time of emplacement, cooling rates of different stages, exhumation/erosion rates for country and intrusion rocks, eroded thicknesses of intrusions, and the exposure time if the sample was exhumed to the surface (McInnes et al., 2005a;Fu et al., 2010).

Zircon U-Pb ages
Most of the zircon grains from the Pulang complex are colorless and transparent, euhedral with clear concentric and oscillatory growth zoning in cathodoluminescence (CL) images. No inherited zircon or relict cores were observed in this study (Electronic App. Table A1). Three indistinguishable U-Pb ages (217.6 ± 1.3-216.4 ± 1.4 Ma) were obtained for the quartz monzonite porphyry samples collected at different sites of the complex (Table 2; Fig. 3a-c). Five quartz diorite porphyry samples also yielded similar zircon U-Pb ages, ranging from 218.1 ± 1.6 to 213.5 ± 1.9 Ma (Fig. 3d-h). These ages represent the maximum emplacement timing of the complex.
(U-Th)/He ages In general, three to five single zircon and apatite grains were analyzed for each sample. The detailed results are given in Electronic Appendix Tables A3 and A4, and a summary with the weighted mean age for each sample is provided in Table  2. Overall, 32 single-grain zircon (U-Th)/He ages, varying from 167.1 ± 9 to 86.7 ± 4.7 Ma, were obtained from the samples. Notably, the age variations between single crystals from some samples (e.g., PLZK2812-1196.1m, PLZKE002-3.5m) exceed the analytical precision of <6%. In addition, 36 single-grain apatite (U-Th)/He ages, ranging from 156 ± 16.2 to 24.3 ± 0.9 Ma, were obtained. Similarly, there are also large intrasample age variations, which are beyond the analytical uncertainty of <2.5%. In order to avoid the potential influence of outliers on averaging dates of each sample, the Chauvenet criterion was utilized. This method has been used for outlier detection in previous studies (e.g., Fitzgerald et al., 2006;Liu et al., 2017). Consequently, seven single-grain zircon (U-Th)/He ages and nine single-grain apatite (U-Th)/He ages that failed the Chauvenet's criterion test were excluded from the determination of the mean and weighted mean ages (Electronic App. Tables A3, A4). Although the single-grain ages are relatively scattered, the weighted mean zircon and apatite (U-Th)/He ages generally display a positive correlation with elevation (Table 2; Fig. 4), allowing us to estimate long-term erosion rates of the Pulang complex. Using a least squares regression routine, we obtained two regression lines to fit the age-elevation data. The two lines predict two different exhumation processes, with an apparent erosion rate of 33.6 m/m.y. during the interval of 142 to 110 Ma and a rate of 15.3 m/m.y. from 106 to 54 Ma.

HeFTy modeling
Thermal history modeling results performed by HeFTy v.1.9.3 are plotted in Figure 5. Overall, the results reveal coherent time-temperature paths for all the analyzed samples. At least three phases of cooling were obtained based on the variable gradients of the weighted mean cooling paths: a phase of rapid cooling at a rate of ~80° to 120°C/m.y. in the Late Triassic, a period of slow cooling at a rate of ~3° to 5°C/m.y. from the Late Triassic to Early Cretaceous, and a phase of extremely slow and monotonic cooling at a rate of <1°C/m.y. from the Early Cretaceous to the present day. Additionally, as shown in Figure 5, all the analyzed samples experienced prolonged periods within the helium partial retention zones for both the zircon and apatite (U-Th)/He systems. Specifically, the residence time in the apatite helium partial retention zones was greater than 20 m.y., whereas the stay in the zircon helium partial retention zones lasted for ca. 10 to 20 m.y.

4DTherm modeling
Based on geologic mapping and drill core logging at the Pulang deposit, the quartz monzonite porphyry was modeled as a 1.8-by 0.7-km oval cylindrical porphyry stock with a residual height of 3.0 km, while the quartz diorite porphyry was simplified as a 3.7-by 3.7-km cylindrical porphyry stock with a residual height of 4.5 km. The crystallization temperature of the zircon in the Pulang complex was estimated to be ~750°C by Wang (2008). The closure temperatures of zircon and apatite (U-Th)/He systems were calculated by the CLO-SURE program of Brandon et al. (1998) at a cooling rate of 1°C/m.y. (see Electronic App. Tables A3, A4). The initial temperature of felsic magma at Pulang was assumed to be 800°C. Considering that the research area is sited in a tectonically active belt, the basal heat flow is set at 100 MW/M 2 , and the geothermal gradient of country rocks was set as 40°C/km. Other initial parameters are copied from the default values of the 4DTherm v.1.2 packages (see Fu et al., 2010, for details). Four of the analyzed samples were successfully tested (Table  3; Fig. 6). The inferred emplacement depth of the Pulang complex is between 5.0 and 6.5 km, which is broadly consistent with the depth inferred from the Al-in-Hornblende barometer (84-215 MPa, Leng et al., 2018). Average cooling rates for the magmatic and exhumation stages were estimated as 112° to 140°C/m.y., and from 1.18° to 1.4°C/m.y., respectively. The duration of the 500° to 300°C interval was calculated as 0.56 to 1.58 m.y. for the quartz monzonite porphyry, which is consistent with the results obtained using hornblende and biotite 40 Ar/ 39 Ar plateau ages of Wang (2008) (see Table 1).
Moreover, 4DTherm estimated an exhumation/erosion rate of 2.2 to 16.8 m/m.y. for the Pulang complex, which is much lower than that of the country rock (32.7-45 m/m.y.). These rates are roughly commensurate with the values derived from the age-elevation relationships of the apatite and zircon (U-Th)/He data (see above). It is estimated that the Pulang complex was exposed to paleosurface at ~74.3 to 48.9 Ma, and the eroded thicknesses of the complex vary from 558 to 1,099 m.

Evaluation of the thermochronological data
As mentioned above, the variation of (U-Th)/He ages within single grains-especially for apatite (U-Th)/He ages-exceeds relevant analytical error. It is thus essential to investigate potential factors causing the age dispersion. The factors may include variations in grain size, heterogeneous distribution of parental nuclides, implantation of He from neighboring U-and Th-rich minerals, presence of excess 4 He from fluid inclusions and/or undetected U-and Th-bearing particles, radiation damage that changed He retentivity and closure temperature, and cooling rate (e.g., Fitzgerald et al., 2006;Zhao et al., 2015;Danišík et al., 2017;Liu et al., 2017).
If the variable single-grain ages resulted from grain size differences, larger grains should have older ages. However, there is no positive correlation between equivalent spherical radius (R eq) and (U-Th)/He ages, suggesting that grain size is not the dominant factor accounting for the intrasample dispersion in single-grain ages (Fig. 7a, b). The possibility of He implantation can also be ruled out, because no U-or Th-rich minerals were observed near zircon or apatite in thin section. A good linear correlation is seen between corrected He content and effective uranium (eU) concentration for both zircon and apatite (Fig. 7c, d), except for some outliers, arguing against the presence of excess 4 He. Although high radiation dosage of zircon could significantly increase He diffusion rates at low temperature (e.g., Reiners, 2005;Danišík et al., 2017), the relatively low U concentrations (mostly less than 500 ppm) as determined by laser ablation-inductively coupled plasmamass spectrometry (LA-ICP-MS) analyses (Electronic App. Table A1) imply that the accumulated radiation damage is insufficient to cause the wide scatter of zircon (U-Th)/He ages (cf. Reiners, 2005;Zhao et al., 2015). Moreover, neither U nor Th zonation patterns were observed during the LA-ICP-MS analyses, suggesting that the heterogeneous distribution of parental nuclides is also not the dominant factor that resulted in the wide variation of single-grain ages at Pulang. Therefore, the cooling rate seems the most likely factor responsible for dispersion of single-grain ages. As shown in Fitzgerald et al. (2006), if the cooling rate is less than 3°C/m.y., the variation in single-grain ages is much greater than it would be at higher cooling rates. It was proposed that when samples   Ketcham (2005). The green area in each plot denotes the envelope for all acceptable time-temperature paths with a goodness of fit >0.05, whereas the purple areas encompass all good time-temperature paths with a goodness fit of >0.6. The thick black line enveloped in the pink region depicts the weighted mean cooling path for each plot. The regions of helium partial retention zone for apatite (AHePRZ) and helium partial retention zone for zircon (ZHePRZ) were calculated using the CLOSURE program of Brandon et al. (1998)  are slowly cooled or resident in the helium partial retention zones for protracted periods, diffusion can magnify the differences in He concentrations for different-sized crystals (especially for small crystals). Consequently, the age variation between different crystals becomes greater (e.g., Meesters and Dunai, 2002;Fitzgerald et al., 2006). The modeled results from HeFTy and 4DTherm suggest a protracted and extremely slow cooling since the Early Cretaceous (Figs. 5, 6), probably giving rise to a comparatively greater variation in the single-grain apatite (U-Th)/He ages. In contrast, the moderate cooling from the Late Triassic to Early Cretaceous led to a smaller variation in the single-grain zircon (U-Th)/He ages.

Formation and duration of the Pulang porphyry system
Porphyry deposits are produced by dynamically evolving magmatic-hydrothermal systems (e.g., Sillitoe, 2010). Detailed description of parameters can be found in McInnes et al. (2005a) and Fu et al. (2010) 1 Emplacement depth is defined as the distance from the paleosurface to the top of the igneous body 2 The igneous body has cooled when it reaches the same temperature as the surrounding country rock 3 Magmatic cooling refers to the cooling from emplacement to the cooled state as defined above 4 Exhumation cooling refers to the interval from the cooled state through cooling to surface temperature of 10°C 5 Exhumation rate of intrusion refers to the average exhumation rate for the intrusion from the emplacement to the paleosurface 6 Exhumation rate of country rock refers to the average exhumation rate for the country rocks from deep to the paleosurface 7 Age of exposure is defined as the age when the top of the igneous body is exposed to the paleosurface 8 Estimated eroded thickness refers to the eroded thickness of the intrusion since it was exposed to the paleosurface   Understanding the duration of genetic magmatic-hydrothermal processes is fundamental to exploring the mechanisms of porphyry deposit formation (McInnes et al., 2005a, b;Chiaradia et al., 2009Chiaradia et al., , 2013von Quadt et al., 2011;Buret et al., 2016Buret et al., , 2017. The suite of radiometric ages obtained in this study (Table 2), combined with previously published geochronological and thermochronological data for the Pulang deposit (Table 1), places good temporal constraints on the initiation and duration of magmatic-hydrothermal activity at Pulang. The new zircon U-Pb data (217.6 ± 1.3-213.5 ± 1.9 Ma) provide additional constraints on the timing of emplacement of the Pulang complex. Zircon (U-Th)/He ages (Tc = 183°C, at 10°C/m.y.; Reiners et al., 2004) have been used to temporally constrain the termination of a magmatic-hydrothermal event (e.g., McInnes et al., 2005a, b;Fu et al., 2010;Braxton et al., 2012) and can be combined with zircon U-Pb ages to reveal the duration of hypogene ore formation. However, this approach is not suitable for Pulang. As shown in Table 2 and Figure 4, the zircon (U-Th)/He ages display a wide range  and generally show a positive correlation with elevation, suggesting that the Pulang complex was emplaced at depths below 4.5 to 6.0 km (assuming a geothermal gradient of 30°-40°C/km), where the ambient temperature exceeds the helium blocking temperature for zircon. Our inverse thermal modeling also yielded a similar emplacement depth (5.0-6.5 km) for the complex. In this scenario, the zircon (U-Th)/ He age records the cooling time of when the sample was uplifted and cooled through the He closure temperature. The biotite 40 Ar/ 39 Ar system has a Tc of 310 ± 20°C (Harrison et al., 1985), which is similar to the lower-temperature limit of magmatic-hydrothermal activity. We therefore propose that it presents a plausible proxy for bracketing the duration of the magmatic-hydrothermal system for deep-seated porphyry intrusions. Three biotite 40 Ar/ 39 Ar plateau ages from 214.6 ± 0.9 to 210.3 ± 0.8 Ma have been presented by previous workers (i.e., Zeng et al., 2006;Wang, 2008;Li et al., 2011; Table  1). Combined with the new zircon U-Pb ages, it is estimated that magmatic-hydrothermal activity lasted at least 3 to 7 m.y. This long time scale is also supported by inverse thermal modeling results (Table 3;  that deposits may form within tens of thousands of years (e.g., Henry et al., 1997;Shinohara and Hedenquist, 1997;Pollard et al., 2005;von Quadt et al., 2011;Chiaradia et al, 2013;Buret et al., 2016Buret et al., , 2017 or several million years (e.g., Ballard et al., 2001;Maksaev et al., 2006;Harris et al., 2008;Sillitoe and Mortensen, 2010;Barra et al., 2013;Deckart et al., 2013, 2014, andreferences therein;Leng et al., 2013). These longlived magmatic-hydrothermal events are typically attributed to episodic buildup and evolution of the underlying more voluminous pluton that feeds the upper crustal porphyritic intrusive complex (e.g., Sillitoe, 2010;Sillitoe and Mortensen, 2010;Chelle-Michou et al., 2014;Chiaradia and Caricchi, 2017;Leng et al., 2018).

Thermal and erosion history of the Pulang deposit
In this study, two different inverse thermal modeling programs (i.e., HeFTy and 4DTherm) were independently utilized to reconstruct the cooling and exhumation histories of the analyzed samples. In general, the time-temperature cooling paths determined by both programs are consistent (Figs. 5, 6), despite minor differences in inflection points. Moreover, the cooling curves in Figure 6 also match the previous published radiometric ages within error, indicating that the modeled cooling paths are reliable. As a whole, the curves in Figure 6 show a typical hockey-stick pattern that is generally attributed to a slowly cooled, deep pluton (e.g., McInnes et al., 2005a). Our modeling yielded an emplacement depth below 5.0 km. In addition, at least three distinct cooling stages can be recognized from the time-temperature cooling curves (Figs. 5, 6). The first stage of rapid cooling in the Late Triassic was related to cooling of the magma by advective heat loss and convection (Cathles, 1981). This process lasted approximately 3 to 7 m.y. as constrained by the zircon U-Pb ages and modeling results (Figs. 5, 6; Table 3).
Both the second and the third stages of cooling at Pulang were primarily controlled by erosion and exhumation processes. Moreover, the decrease in the cooling rate (i.e., from ~3°-5°C/m.y. in the second stage to less than 1°C/m.y. in the third stage) could be related to a decrease in the erosion rate during the third phase, consistent with the results derived from the age-elevation relationships (Fig. 4). Collectively, these results imply that Pulang has undergone at least two stages of exhumation (Fig. 8). The early-stage erosion was initiated no later than 182 m.y. ago based on the K-feldspar 40 Ar/ 39 Ar plateau age (Wang, 2008), probably as early as the Late Triassic, and lasted until the Late Cretaceous. During this stage, the country rocks were exhumed at a rate of 33 to 45 m/m.y. The estimated thickness of eroded country rock was more than 5.0 km. The Pulang complex was probably exposed to paleosurface between the Late Cretaceous and Early Eocene, based on the results of our modeling. The late stage continued to the present day, unroofing the Pulang complex at a much slower rate of 5 to 17 m/m.y., with ~558-to 1,099-m thicknesses of materials removed likely. A substantial amount of ore from the Pulang deposit was probably removed during uplift and erosion, given that the Cu-bearing stockwork veins are exposed at surface (Fig. 2b). The long-term (>50 m.y.) exposure of the Pulang mineralization, combined with the moderate to strong rainfall in this region, inhibited supergene enrichment of the deposit.

Implications for the preservation of ancient porphyry systems in southeastern Tibet
It is suggested that the entire Yidun terrane, as well as broader southeastern Tibet, cooled below 300°C since the Late Triassic (~199 Ma) and experienced a long period (>100 m.y.) of thermal stability from the Jurassic to early Paleogene (e.g., Reid et al., 2005;Roger et al., 2011;Tian et al., 2014). Our dating results from the Pulang complex provide more detailed insights into the Cenozoic thermal evolution of this region, which involved an early stage of relatively rapid erosion and a late stage of much slower erosion. This decrease in cooling rate in the Yidun terrane was previously considered to be a result of crustal refrigeration due to the subduction of the Meso-Tethanys and the subsequent collision between the Lhasa and Qiangtang terranes in the Cretaceous (e.g., Tian et al., 2014). Kapp et al. (2005) has demonstrated that the Cretaceous collision of the Lhasa and Qiangtang terranes resulted in substantial shortening and thickening of the upper crustal rocks, causing regional uplift of the eastern Tibetan Plateau during the early Paleogene (~45 Ma; Rohrmann et al., 2012). However, this is not the case in the Yidun terrane. Considering the geomorphic features and thermal evolution of the Yidun terrane, we therefore suggest that moderate cooling in the upper crustal rocks from the Late Triassic to Early Cretaceous as recorded by the Pulang complex resulted from isostatic adjustment of the crust after the subduction of the Ganz-Litang Ocean beneath the Yidun terrane in the Late Triassic . A low-relief paleosurface landscape (i.e., peneplain), probably at low altitude, formed during the Late Cretaceous. The slow crustal cooling rate and extremely slow erosion rate from Late Cretaceous to the present day was a response to the presence of this peneplain landscape. A number of low-temperature thermochronological studies from several large rivers (e.g., Yalong, Dadu, and Jinsha Rivers) that dissect the Yidun terrane indicated that rapid river incision of 200 to 500 m/m.y. started around 15 to 9 Ma (e.g., Clark et al., 2005;Ouimet et al., 2010). This rapid incision was considered to be a response to the nearly coeval rapid uplift of the entire southeastern Tibetan Plateau. Clark et al. (2005) invoked a lower crustal and mantle channel flow model to explain the uplift of the southeastern Tibetan Plateau. However, this model was recently challenged by several studies from adjacent areas that demonstrate that the onset of surface uplift commenced from late Oligocene to early Miocene (30-15 Ma; e.g., Wang et al., 2012;Tian et al., 2014). These obvious discrepancies in the uplift history from different sites strongly support a stepwise uplift geodynamic model (Tapponnier et al., 2001). The prolonged period of extremely slow cooling and denudation of the Pulang complex indicates that the porphyry system was not affected by rapid Miocene denudation. Therefore, the Miocene rapid uplift was not propagated to the entire Tibetan Plateau but was likely restricted to some major faults (e.g., Roger et al., 2011;Zhang et al., 2015Zhang et al., , 2017. In summary, we propose that the presence of a relict peneplain in modern southeast Tibetan Plateau was favorable for preservation of the Triassic Pulang porphyry deposit. Our results may also provide new insights into the preservation of other ancient porphyry deposits in the pre-Mesozoic orogenic belts, such as the Devonian Oyu Tolgoi porphyry Cu-Au deposits in the Central Asian orogenic belt (e.g., Wainwright et al., 2017) and the Silurian Cadia porphyry Au-Cu clusters in the Lachlan fold belt (e.g., Holliday et al., 2002).

Concluding Remarks
U-Pb-He dating and numerical modeling suggest that the Pulang intrusive complex was emplaced at a paleodepth of 5.0 to 6.5 km around 215 ± 2 Ma. The relatively deep level emplacement of the complex, coupled with the episodic replenishment of the magma chamber, accounts for prolonged magmatic-hydrothermal activity at Pulang.
Three phases of cooling were revealed by the thermal history modeling at Pulang, including a rapid cooling phase in the Late Triassic, a moderate cooling from the Late Triassic to Early Cretaceous, and a monotonic and extremely slow cooling from the Early Cretaceous to the present day. The early, first phase of cooling mainly reflects magmatic cooling processes, whereas the later two phases of cooling were predominantly controlled by erosion and exhumation processes. Moreover, the decrease in cooling rate from the second to the third phase could be related to a decreasing erosion rate. Overall, our results indicate that the Pulang complex has experienced at least two stages of exhumation at rates of 33 to 45 m/m.y. and 5 to 17 m/m.y., respectively. Moreover, it is estimated that a substantial amount of ore from the Pulang deposit was probably removed during uplift and erosion processes.
The Triassic Pulang deposit has undergone long-term (>50 m.y.), extremely slow cooling (<1°C/m.y.) and denudation since the Late Cretaceous and was not affected by the rapid incision in the late Miocene. This is probably related to the formation of a relict peneplain on the southeastern Tibetan Plateau from the Cretaceous that has survived to the present day.