We performed apatite and zircon (U-Th)/He dating on a granitic pluton that has been offset by ∼10 km by motion on the sinistral strike-slip Xiangcheng fault in SW Sichuan, SE Tibetan plateau, where the Shuoqu River incises a deep valley before joining the upper Yangtze River. Mean ZHe cooling ages range from 49.5 ± 2.2 Ma to 68.6 ± 6.0 Ma. Samples located above 3870 m yield mean apatite (U-Th)/He ages ranging from 30.6 ± 1.4 Ma to 40.6 ± 2.7 Ma, whereas samples at lower elevations range from 9.8 ± 1.3 Ma to 14.6 ± 2.7 Ma. In the same region, Cenozoic continental sediments are exposed on the flanks of deep valleys. They consist of unsorted conglomerates and sandstones that partly fill a paleotopography. The sediments were deposited during an episode of rapid sedimentation, followed by incision that varies between 0.5 and 1.2 km. Thermal and exhumational modeling of the granite thermochronometric data indicates rapid cooling during the middle Miocene that was likely related to fluvial incision. Our findings suggest that the upper Yangtze River and its tributary (Shuoqu) were connected by the middle Miocene. Our modeling also supports the idea that the exhumation pattern during the Cenozoic in the southeastern margin of the Tibetan Plateau is spatially and temporally heterogeneous.

The landscape on the southeastern margin of Tibet is characterized by narrow and deeply incised fluvial valleys (Fielding et al., 1994) with sharp bends close to faults. It transitions northward to the low-relief areas of the Tibetan Plateau and southwards into the lower-elevation reaches of Yunnan (Liu-Zeng et al., 2008). This landscape formed during the India-Asia collision in response to tectonic deformation and uplift (Brookfield, 1998; Clark and Royden, 2000). Several geodynamic models propose different mechanisms to explain the plateau formation, including homogeneous thickening (England and Houseman, 1986), lateral extrusion (Tapponnier et al., 2001), and lower-crustal flow (Royden et al., 1997). The landscape response to collision has been complex, including constriction of drainage basin geometry (Hallet and Molnar, 2001), major river pattern reorganization (Clark et al., 2004; Yang et al., 2015), and progressive incision possibly reflecting late uplift (Clark et al., 2006). A key to understanding the processes associated with the formation of this landscape is establishing the timing and amount of fluvial erosion. The southeastern margin of Tibet is drained by three main rivers: the Salween (Nu Jiang), Mekong (Lancang Jiang), and Yangtze (Chang Jiang; Fig. 1). Several studies based on different techniques have constrained erosion in several locations along these rivers at the million-year time scale (Clark et al., 2005; Ouimet et al., 2010; Wilson and Fowler, 2011; Duvall et al., 2012; Tian et al., 2014; Ping et al., 2015; Zhang et al., 2015; Shen et al., 2016; Yang et al., 2016; McPhillips et al., 2016; Liu-Zeng et al., 2018; Nie et al., 2018) and at the 1–100 k.y. time scale (Henck et al., 2011; Zhang et al., 2018). However, no clear pattern or single mechanism has emerged (see review by Zhang et al., 2019).

The Yangtze River is particularly intriguing because it has a complex course (Figs. 1 and 2A): Along its upper portion (Jinsha River), it flows southward and parallel to the Mekong and Salween Rivers, and then it bends abruptly in the Yunnan highlands. There, it forms a V-shaped loop, called the Yangtze first bend, and becomes the middle Yangtze River flowing northward. It then bends again in northern Yunnan at latitude ∼27.75°N to flow southward and then changes once again to flow eastward and northeastward along the SE edge of the Sichuan Basin. Along its last portion, the lower Yangtze flows to the East China Sea. The paleocourse of the Yangtze River has been discussed for more than a century (e.g., Abendanon, 1908; Gregory, 1929; Ting, 1933; Lee, 1934; Barbour, 1936; Ren et al., 1959; see review by Zheng, 2015). Recent geology studies on that topic are often based on thermochronometry and detrital sedimentology and aim at constraining the timing of fluvial incision and establishment of the river course.

Based on isotopic detrital data along the middle Yangtze River, Zhang et al. (2017) placed a late Pliocene upper boundary for the establishment of the modern river, whereas Zheng et al. (2013) suggested it occurred before the Miocene, based on detrital zircon age data along the lower Yangtze and near the Yangtze outlet. One of the main questions regarding the evolution of the Yangtze is a potential former connection between the upper Yangtze and the Red River, which would correspond to a paleo-Yangtze flowing southward (e.g., Clift et al., 2006; Yan et al., 2012; Wei et al., 2016; Wissink et al., 2016; Chen et al., 2017; Gourbet et al., 2017). Near the Yangtze first bend, the Jianchuan Basin (JB in Fig. 2) contains fluvial sediments of Eocene age (Gourbet et al., 2017), which are interpreted as relicts of the Yangtze River (Clark et al., 2004). Provenance reconstructions based on detrital zircon ages from these sediments and from sediments near the Red River upper reaches indicate a main river flowing from the plateau interior southward during the Paleogene (Yan et al., 2012; Chen et al., 2017). These reconstructions are consistent with the hypothesis of major river modifications at 35–31 Ma based on isotopic data from paleo–Red River sediments in the Hanoi Basin (Clift et al., 2006). Other reconstructions, based on detrital zircon ages and on paleoflow directions from several basins, including the Jianchuan, Lijiang, and Yanyuan basins (Fig. 2), indicate no Cenozoic connection between the paleo–Yangtze River and paleo–Red River (Wei et al., 2016; Wissink et al., 2016). The truncation of a hypothetic southward-flowing paleo–Yangtze River might have been due to the activation of the Ailao Shan–Red River fault during the late Eocene (Gourbet et al., 2017; Chen et al., 2017), thus leading to the modern upper Yangtze and Red River configuration. Other potential triggers for drainage reorganization are onset of regional surface uplift (Clark et al., 2004) or local volcanic activity (Zeng, 2002; Gourbet et al., 2017). Today, portions of the course of the upper Yangtze follow active strike-slip faults. The river is offset by ∼80 km by the left-lateral Xianshuihe fault (XF in Fig. 1; Gaudemer et al., 1989). Further south, the Batang fault right laterally deflects the river course for ∼30 km (Fig. 2; Wang and Burchfiel, 2000). The Yangtze is also offset with a left-lateral sense for ∼15 km by the Zhongdian fault (Fig. 2; Burchfiel and Wang, 2003), although this fault may have experienced a reversal in sense of slip (Perrineau, 2010).

Low-temperature thermochronometry can address the timing of incision and thus potential reorganization events. However, many terranes along the main trunk of the Yangtze have lithologies with only rare apatite and zircon, which limits the number of sites suitable for thermochronometric studies. Along the main upper Yangtze valley, Wilson and Fowler (2011) reported five apatite fission-track ages ranging from ca. 12 to 6 Ma, and Ouimet et al. (2010) reported four apatite (U-Th)/He (AHe) ages from an elevation transect ranging from ca. 9 to 7 Ma. In tributaries located a few kilometers west of the first bend, moderate river incision occurred from 18 to 9 Ma (McPhillips et al., 2016). However, in the same area, Yang et al. (2016) reported a late Eocene AHe age. Based on AHe ages ranging from 30 to 20 Ma, Shen et al. (2016) proposed that incision occurred during the Miocene–early Oligocene, which was confirmed by new AHe and fission-track data by Cao et al. (2019), who showed a rapid exhumation phase from 28 to 20 Ma. Downstream of the Yangtze first bend, in the Tiger Leap Gorge (TLG in Fig. 2), Lacassin et al. (1996) suggested that incision by the Yangtze was ongoing around 17.5 Ma during folding of an anticline, based on K-feldspar Ar-Ar thermochronometry.

Here, we focused on one of the tributaries of the upper Yangtze River, the Shuoqu River. It flows southward from the SE Tibetan Plateau into the upper Yangtze River (Fig. 2). About 70 km east of the Yangtze and over a distance of 40 km, the Shuoqu River flows along a granitic pluton and Cenozoic continental sediments, crossing a strike-slip fault (called the Xiangcheng fault). In this study, we applied apatite and zircon (U-Th)/He (ZHe) thermochronometry to estimate cooling rates, which we used to quantify the erosion in the area. We combined these data with new observations of perched detrital sediments. Ultimately, we compared our results with published studies and extended our discussion to the regional geological evolution of the southeastern Tibetan Plateau.

The Xiangcheng area is located in SW Sichuan, at the transition between the high-elevation, low-relief landscapes of the eastern Tibetan Plateau and the Yunnan highlands (Figs. 1A and 1B). Rivers flow either approximately SW toward the upper Yangtze reach, as is the case for the Shuoqu River, or approximately S/SE toward the middle Yangtze reach (Figs. 2 and 3). These rivers are more than 100 km long, in contrast to the short (10–20 km) eastward-flowing tributaries of the upper Yangtze reach, which drain the narrow area confined between the Mekong and Yangtze Rivers. Local relief around the city of Xiangcheng has a magnitude of nearly 2.5 km and a wavelength of >20 km (Fig. 2). In this area, the Shuoqu River is a dominantly bedrock river associated with a few quaternary alluvial terraces (Fig. 4).

The study area belongs to the Yidun arc (or Litang-Batang block; Fig. 2). To the west, it is separated from the Qiangtang block by the Jinsha suture zone (Şengör, 1985). To the east, it is separated from the Songpan-Garzê block by the Garzê-Litang suture zone. The western Yidun arc is mostly composed of Paleozoic carbonates and volcanic rocks (Pan et al., 2004). It is considered to be a remnant of a passive margin that was formerly attached to the northern South China block. The eastern Yidun arc, where the study area is located, is composed of Triassic continental and marine sandstones, slates, and limestones. Magmatic rocks are Late Triassic calc-alkaline volcanic-to-subvolcanic and plutonic rocks, and Cretaceous granites. Cenozoic continental sediments are rare, but they outcrop in the Yidun arc (Fig. 2; Pan et al., 2004). We describe these sediments in the “Sedimentology and Paleolandscape” section.

The Yidun arc experienced a complex deformation history. The closure of the Paleotethys ocean due to convergence of the Qiangtang, North, and South China blocks led to three subduction zones, including the Litang-Batang subduction (to the east) and the Jinsha subduction (to the west), which were active during the Early-Middle Triassic (Wang et al., 2000). The study area is now sandwiched between the Qiangtang and the Songpan-Garzê blocks (Fig. 2). It is affected by approximately N-S–striking folding and thrusting related to the formation of an accretionary Triassic orogenic belt (Reid et al., 2005; Roger et al., 2008). Near the city of Xiangcheng, folds have an approximately SSE-NNW axis and a kilometer-scale wavelength (Reid et al., 2005). The end of this deformation stage is characterized by undeformed Late Triassic and Cretaceous plutons intruding the deformed sediments (Reid et al., 2005, 2007), including a granitic pluton we sampled for thermochronometry dating near the city of Xiangcheng (see “Sampling and Methods” section). This granite belongs to the so-called “Cilincuo intrusion” (Wang et al., 2008; Fei et al., 2009; Li et al., 2014; Wang et al., 2014), for which previous studies reported U-Pb ages ranging from the Late Jurassic to the Paleocene. Cu deposits in this intrusion were dated at ca. 79 Ma (U/Pb; Fei et al., 2015), whereas porphyritic granites with amphibole and biotite were dated at 147 ± 4.6 Ma (U/Pb on zircon; Wang et al., 2008). However, Wang et al. (2008) did not provide the geographical coordinates of the samples, and the location of this early intrusion is uncertain. Last, Li et al. (2014) quoted ages of 98 and 51 Ma from the Yunnan Geological Survey without providing the dating method nor the sample locations. Thus, the Cilincuo intrusion could actually consist of multiple intrusions, and in order to avoid ambiguities, we define the Cilincuo massif as the granitic massif comprising these intrusions.

After a period of relatively limited tectonic activity during the Jurassic–Cretaceous (Roger et al., 2010), deformation resumed following the India-Asia collision and subsequent E-W shortening. Cenozoic deformation in SE Tibet is characterized by EW/SW-NE shortening (Lacassin et al., 1996) and active approximately NW-SE strike-slip faulting (Tapponnier and Molnar, 1977). The Cenozoic shortening affected continental basins, which show contractional structures coeval with sedimentation (Horton et al., 2002), and possibly consisted of two distinct phases at ca. 36 Ma and ca. 17 Ma (Lacassin et al., 1996). In the Mula Basin (MB in Fig. 2A), Jackson et al. (2018) reported NE-SW–oriented thrust faults that are either at the basin boundary or in the basin and affect sediments that are younger than ca. 45.5 Ma and possibly early Oligocene in age. In the late Eocene Jianchuan Basin (JB in Fig. 2A), sediments are affected by folds with an approximately N-S axis consistent with approximately E-W shortening (Bureau of Geology and Mineral Resources [BGMR], 1991; Gourbet et al., 2017).

The main active strike-slip faults in the area are the left-lateral Litang fault in the Yidun arc (Wang and Burchfiel, 2000; Zhang et al., 2015) and the major Red River fault further south (Leloup et al., 1993; see Fig. 2). The latter was involved in the lateral extrusion of continental blocks (first the Indochina block and then the South China block; Tapponnier et al., 2001). The left-lateral Xiangcheng fault extends along a distance of ∼50 km (Fig. 3). It is not clear where it ends, especially toward the southeast portion. The Shuoqu and Mayi Rivers seem to be left-laterally deflected along the fault, which may indicate that the fault is active. Moreover, two recent ∼10-km-deep, 3–4 Ms earthquakes (Fig. 3) have occurred in the vicinity of the fault (csndmc.ac.cn; earthquake.usgs.gov). However, as microseismicity is ubiquitous in the Tibetan Plateau, these earthquakes do not prove the fault activity, which remains unclear.

About 40 km south of Xiangcheng, Cenozoic clastic sediments, or red beds (Pan et al., 2004), are preserved near the village of DongWang (site 4, Fig. 3). The term “red beds” is widely used in the literature to describe Cenozoic sediments in Tibet; here, we only use this term to refer to the ones in the Xiangcheng area. The red beds are exposed on the flanks of deeply incised valleys, where most of the bedrock is composed of intensively folded Triassic marls and limestones (Fig. 3). The Cenozoic red beds consist of medium-grained sandstones and conglomerates. From bottom to top, the strata show an alternation of red conglomerates (Ed1), sandstones (Ed2), and conglomerates (Ed3). Sediments are slightly deformed; they are gently folded, but in some places, they are mapped as horizontal (BGMR, 1991). We observed planar sandstone bedding in places. Thrust-related deformation affects the sediments over a few tens of meters (Fig. 5A). The contact between the red beds and the underlying Triassic bedrock corresponds to an angular unconformity, above which the red beds unconformably rest on Triassic marls and partly fill a paleorelief that existed prior to deposition of the red beds (Fig. 5C). The sedimentary unconformity appears at several locations in the DongWang area (Fig. 4). However, according to BGMR (1991), the westernmost red beds are bounded by a west-dipping, N-S–oriented thrust fault, with the Triassic units overthrusting the Cenozoic red sediments (Fig. 4). Tectonic control on Cenozoic sedimentation has been reported elsewhere in SE Tibet, including the Mula Basin (Jackson et al., 2018). Thus, it is possible that the deposition of the red beds was (at least partly) controlled by the local thrust fault identified by BGMR (1991). Today, the red beds are cut by the Shuoqu River and its tributaries. Incision is highly variable if it is inferred from the elevation difference between the top unconformity of the Cenozoic deposits and the valley bottom. For example, directly above the hamlet in Figure 3, the incision is ∼500 m, whereas it can be up to 1200 m further east.

The red beds belong to the Dongwang formation, which is attributed to the Eocene–Oligocene (Zhang, 2009; Yunnan Regional Geological Survey Team, 1982). According to BGMR (1991), minimum thicknesses of Ed1, Ed2, and Ed3 sediments are 237, 288, and 479 m, respectively. Conglomerates consist of a red sand-gravel mixture matrix, unsorted clasts, and boulders of a few centimeters to decimeters in size. Boulders are angular to subrounded, suggesting limited transport, and are apparently solely composed of limestone (Fig. 5D), likely originating from the surrounding Triassic peaks or at least coming from a local source. These sediments do not correspond to a large riverbed, but rather indicate a short episode of intense sedimentation with a significant contribution from the talus, in a fluvial environment. The cause of this episode is unclear and may have been related to the activation of the local thrust fault mentioned above. This hypothesis is consistent with a local sedimentary source. The Triassic marls and limestones are well constrained in age because of their fossil record (BGMR, 1991). However, to our knowledge, there is no absolute nor relative age constraint on the Cenozoic red beds that would allow us to directly infer an incision rate from the incision amount we observed in these sediments. Because Cenozoic continental sediments are difficult to date, the difference between the assumed and revised depositional age in SE Tibet can reach more than 20 m.y. (Gourbet et al., 2017; Linnemann et al., 2018). As a consequence, we cannot estimate an incision rate based on the age of the red beds, and so we simply used the observed range of 0.5–1.2 km as an estimate for the total incision since deposition.

Sampling Strategy

We collected granite samples mainly at three locations in the Cilincuo massif (Figs. 2 and 4; Data Repository Table ST11): site 1 at its northern edge, site 2 directly south of the strike-slip Xiangcheng fault, and site 3 at its southern edge. For site 1, samples were distributed over a distance of 1.7 km and ∼160 m of elevation; for site 2, samples were distributed over a distance of 2.6 km and ∼270 m of elevation; for site 3, samples were distributed over a distance of 2 km and ∼550 m of elevation. We also collected a sandstone sample (S16X14) located near the valley floor along the main N-S road, ∼8 km west of the intrusive contact, for thermochronometric dating. We collected an additional single granitic sample, S16X19, further south in site 4 for U-Pb and (U-Th)/He dating.

U-Pb Dating of Zircon

U/Pb dating was performed by laser ablation–inductively coupled plasma–mass spectrometry (ICP-MS) at the Laboratoire Magmas et Volcans, Clermont-Ferrand, France. The analytical results are presented in the Data Repository (Table ST1). Zircons were ablated using a Resonetics Resolution M-50 laser system operating at a wavelength of 193 nm coupled to a Thermo Element XR ICP-MS. Helium carrier gas was supplemented with N2 prior to mixing with Ar for sensitivity enhancement (Paquette et al., 2014). The laser was operated with a spot diameter of 27 μm, a repetition rate of 3 Hz, and a fluence of 2.5 J/cm2. Instrumental operating conditions and data acquisition parameters are basically similar to those reported in Hurai et al. (2010). Reduction of raw data was carried out using the GLITTER® software package (Van Achterbergh et al., 2001). Isotope ratios were corrected for laser-induced and instrumental mass fractionation via sample-standard bracketing using the GJ-1 zircon (Jackson et al., 2004). Data were not corrected for common Pb. The 207Pb/206Pb versus 238U/206Pb diagram was generated using the Isoplot/Ex v. 2.49 software of Ludwig (2001). Error ellipses for each point are shown at the 2σ level and incorporate both internal and external uncertainties. Data points were pooled to calculate a date and associated 2σ error. The 91500 zircon reference material (Wiedenbeck et al., 1995) was analyzed along with the samples to independently monitor the external precision and accuracy of the measurements. The concordia age for 132 analyses of 91500 zircon conducted over the course of the study was 1063.9 ± 2.4 Ma (2σ, including decay constant errors).

(U-Th)/He Thermochronometry

The (U-Th)/He method allows assessment of the timing of cooling of a rock from relatively low temperatures to ambient surface conditions. For apatite, the nominal closure temperature (Tc) of the (U-Th)/He system is ∼65 °C (Shuster et al., 2006). The partial retention zone (PRZ) for apatite spans from 40 °C to 85 °C (Wolf et al., 1998). He diffusivity can be altered by several factors, including the accumulation of radiation damage, which in turn causes Tc to span from 45 °C to 120 °C (Shuster et al., 2006; Flowers et al., 2009). For zircon, Tc is ∼180 °C at a cooling rate of 10 °C/m.y. (Reiners, 2005) and varies from 140 °C to 220 °C depending on the self-irradiation dose (Guenthner et al., 2013); the PRZ for zircon spans from 80 °C to 180 °C (Ehlers et al., 2005). Thus, combining apatite (U-Th)/He and zircon (U-Th)/He (AHe and ZHe) dating provides information on the exhumation processes in the upper 6 km of the crust at medium geothermal gradients (∼30°/km).

All (U-Th)/He analyses were performed at the Institute of Geochemistry and Petrology, ETH Zürich, and are summarized in Table 1. We attempted to obtain four ages from four grains for each sample. In three cases, we were unable to find sufficient apatite of high quality (Table 1). A complete description of sample preparation and analytical technique is provided in Yang et al. (2016), and the error and reproducibility calculations are explained in the Data Repository.

Zircon U-Th Zonation

Hourigan et al. (2005) demonstrated that intracrystalline zonation of U and Th in zircon can lead to an inaccurate estimate of the alpha-ejection correction and ultimately to age inaccuracy. In order to evaluate the potential effects of U and Th zonation on intrasample ZHe age variations, we performed additional U and Th measurements on sample S16X10 by laser-ablation ICP-MS. We used 32 of the 50 grains that were used for U-Pb dating and performed an additional measurement on each crystal rim, and then we compared the U and Th concentrations in the rim with the concentrations already measured in the center of the grains during U-Pb dating. The analytical results are presented in the Data Repository (Table ST2; Fig. S2). We followed the same method as described earlier herein, except that the laser was operated with a spot diameter of 7 μm, thus allowing us to access the most external part of the crystal rim. We present the results in the “U-Pb Zircon Geochronology” subsection of the “Results” section.

Modeling of Cooling and Erosion

We inferred the thermal history of the Cilincuo massif using a one-dimensional (1-D) inverse modeling approach, in order to derive thermal history information from thermochronometric data. We used the published code, QTQt (Gallagher, 2012), which uses thermochronometric kinetic models to find an optimum temperature history for multiple chronometers. It also permits the simultaneous inversion of multiple samples, provided they have a common exhumation history but differ in depth, and so have a constant temperature offset in time, as occurs with near-vertical sample profiles. The list of input and output data is provided in Table 2. For each sampling site, we modeled granite samples together: samples S16X15, S16X16, and S16X17 for site 1, samples S16X01, S16X03, S16X04, S16X05, and S16X08 for site 2, and samples S16X09, S16X09B, S16X10, S16X11, S16X12, and S16X13 for site 3. Because several samples from the same site but with different elevations were modeled simultaneously, each run required input of the temperature offset corresponding to the temperature between the sample with the highest elevation and the sample with the lowest elevation, i.e., depending on the geothermal gradient. For site 1, the elevation range between the topmost and lowermost sample was 163 m; for site 2, it was 268 m; for site 3, it was 552 m. For the final temperature offset (i.e., for a given site, the temperature difference between the topmost and lowermost sample once they are at the surface), we used the surface temperature difference calculated using a 5–6°/km adiabatic gradient (Table 2). Thus, for each simulation, the temperature offset varied from the value corresponding to the initial geothermal gradient to the surface temperature offset.

The program QTQt samples the temperature-time (T-t) space into numerous temperature-time (Ti, ti) points, in order to discretize the time-temperature space. Then, using a Markov chain Monte Carlo method, it tests multiple cooling models to find the best cooling history that (1) fits modeled ZHe and AHe ages and (2) limits the total number of T-t points in the history, thereby avoiding overly complex thermal histories. For each T-t path, single ZHe and AHe ages are predicted using helium diffusion kinetics from Ketcham et al. (2011). The fit between the observed ages and predicted ages is calculated using the likelihood function (Gallagher, 2012):
with N being the number of iterations and σ being the observed age error. Then, a decision is made either to retain the simulation if the fit is high or discard it if it is low. This process is repeated until an adequate set of solutions is obtained or the best solution is regarded as adequate. We tested different parameters, including the number of iterations and number of T-t points to reach a thermal history we considered as acceptable, or until further iterations failed to provide a significantly improved model. We found that 20,000–40,000 iterations were sufficient to find convergent behavior. Second, we tested the influence of the geothermal gradient on the thermal histories. In the “Thermal Modeling” subsection of the “Interpretation and Discussion” section, we describe our preferred models, whereas in Data Repository Section 5, we present other models with different input parameters.

Field Observations

Lithology in the Cilincuo massif is variable. The main phase includes quartz, biotite, feldspar, and amphibole. The granite texture is often porphyritic (Wang et al., 2008) with centimeter-scale pink/white K-feldspar crystals (Data Repository Fig. S1). We observed the intrusive contact between the granite and Triassic units at site 1 (Data Repository Fig. S1). Sampled outcrops are fresh, but in other places, the granite is highly weathered and arenitic. In site 2, the granite is cut by an ∼10-cm-wide aplitic dike that we sampled for both U/Pb and (U-Th)/He dating (sample S16X01). The granite intrudes folded Triassic sediments that consist of pelites from the Lamaya formation (BGMR, 1991). Lower in the sedimentary succession, sediments consist of interbedded sandstones, marls, and pelites (BGMR, 1991). At site 4 (sample S16X19), the granite is weathered, and thin sections show chloritized feldspar and calcite recrystallization. The apatites from this sample, observed under the microscope after mineral separation, are weathered.

We mapped the Xiangcheng fault using the Cilincuo massif geometry and the deflection of the courses of the Shuoqu and Mayi Rivers (Fig. 2). In addition, 15-m-resolution LANDSAT images revealed the fault trace along an ∼2-km-long portion of the fault but did not indicate the slip sense. The offsets of the rivers provide minimum values for the total displacement along the fault (Gaudemer et al., 1989). These values are 9 km for the Shuoqu River and 14 km for the Mayi River. The fault also offsets the Cilincuo massif by ∼10 km.

U-Pb Zircon Geochronology

Here, we report results of zircon U-Pb dating of the Cilincuo massif and for granitic sample S16X19, which was collected in another pluton further south at site 4.

We analyzed 37 grains from aplite sample S16X01 (site 2, Cilincuo massif; see Data Repository Table ST2). On the Tera and Wasserburg diagram (Fig. 5A), 36 of the 37 analytical points define a lower intercept at 83.11 ± 0.57 Ma, which corresponds to the crystallization age of the aplite dike and dates the end of plutonism. Fei et al. (2015) described a sliver of monzogranite intruding the granite we sampled for thermochronometry. They obtained a zircon U-Pb age of 79.33 ± 0.7 Ma, which is generally consistent with our data.

From granite sample S16X10 (site 3, southern part of the Cilincuo massif), the 50 analyzed grains exhibit a lower intercept at 83.70 ± 0.47 Ma on the Tera and Wasserburg diagram (Fig. 5B). This age is statistically indistinguishable from that of the aplite dike (sample S16X01) that cuts the massif at site 2. Therefore, when interpreting the thermochronological data set, we used 83 Ma as the crystallization age for all samples.

For granite sample S16X19 (site 4), most of the zircon grains yielded Late Triassic ages, and a linear regression through 20 points gave a lower intercept at 206.9 ± 2.2 Ma (Fig. 5C). This time is consistent with the abundant Late Triassic magmatism in the Yidun arc (Pan et al., 2004). The remaining six grains yielded Neoproterozoic and Late Jurassic ages, which are older than the regional magmatism, and, thus, we interpret them as inherited zircons (concordia diagram in Fig. 5D).

(U-Th)/He Thermochronology of Zircon and Apatite

Mean ZHe ages of the Cilincuo massif are between 68.6 ± 6.0 Ma and 49.5 ± 2.2 Ma (Table 1). The AHe ages from site 2 can be divided into two groups (Fig. 6). The two samples located above 3875 m yielded mean AHe ages of 30.6 ± 2.8 Ma and 40.6 ± 4.7 Ma, and the rest of the samples at lower elevations had ages between 14.6 ± 5.3 Ma and 9.8 ± 2.7 Ma. AHe ages of site 1 ranged from 14.6 ± 5.3 Ma to 11.4 ± 1.1 Ma, and AHe ages of site 3 ranged from 10.3 ± 2.3 Ma to 7.9 ± 2.6 Ma (Table 1).

Triassic sandstone sample S16X14, collected near the valley bottom, yielded a mean Late Jurassic ZHe age of 146.2 ± 8.3 Ma and a mean AHe age of 21.1 ± 6.9 Ma, which is younger than the depositional age of this sample, indicating that the sediments were buried deep enough to reset both apatite and zircon (U-Th)/He thermochronometers.

For most zircon grains of sample S16X10, there was limited difference in Th concentrations between the center and the rim (Data Repository Table ST3; Fig. S2); the ratio Thrim/Thcenter varied from ∼0.5 to ∼3, except for one grain that had a 20-fold–enriched rim. U zonation was also limited, with Urim/Ucenter ranging from 0.5 to 5 for 30 zircons out of 32. The two remaining grains showed a 10-fold U enrichment of the rim.

Cooling Age versus Crystallization Age

The two samples from the Cilincuo Massif dated with the U-Pb method indicate that the crystallization age of the sampled granites is 83 Ma. Out of thirty-eight (38) ZHe ages obtained from these granites, two replicates from samples S16X05 yielded ZHe ages larger than the crystallization age, and we consider these two data points as outliers. The remaining ZHe grain ages from the same granitic rocks varied between 82 and 41 Ma, and such a wide range could reflect several factors affecting He diffusion in zircons, which we discuss in more detail in the next section. Overall, the overlap in the Cilincuo massif between some cooling ages and the crystallization age indicates that it is difficult to discriminate between postintrusion and exhumation-related cooling in the range of temperatures close to the closure temperature of the ZHe system. Exhumation-related cooling is, instead, recorded by the AHe grain ages of the Cilincuo granite, which are between 46 and 8 Ma and are therefore much younger than the crystallization age of the granite.

Sandstone sample S16X14 came from a location ∼5 km distant from the contact with the Cilincuo intrusion, located near the valley floor (∼10 m above the present-day riverbed) at lower elevations than the samples of the Cilincuo massif. The ZHe ages of this sample (176–140 Ma) are much older than the crystallization age and the ZHe cooling ages of the Cilincuo granite intrusion. These observations suggest that the pluton intrusion did not likely affect the sandstone sample and that the sandstone and the Cilincuo granites had different cooling histories at temperatures close to the closure temperature of the ZHe system. It is difficult to say whether, at lower temperatures, the sandstone sample shares a common history with the samples of the Cilincuo massif: In fact, we could only date two apatite grains from this sample, and they gave very different ages (16 and 38 Ma). However, the most important conclusion that we can derive based on this sample is that the rocks near the valley floor have been experiencing temperatures below ∼200 °C for more than 140 m.y. Thus, S16X14 provides an upper limit for the total Cenozoic exhumation along the valley of the Shuoqu River.

For sample S16X19 (site 4), all AHe and ZHe ages are much younger than the ca. 206 Ma crystallization age of this sample. Thus, we interpret the AHe and ZHe data of that sample as exhumation-related cooling. Both the AHe (38–16 Ma) and ZHe (136–68 Ma) cooling ages of this sample are older than or overlap the cooling ages of the samples within the Cilincuo massif, although this sample is located at lower elevations along the main trunk of the Shuoqu River. The very large dispersion of all the cooling ages from this sample does not allow much inference about its cooling history. However, the ZHe ages of this sample are Mesozoic, and its AHe ages are Cenozoic, and, therefore, like sample S16X4, it provides an upper limit for the total Cenozoic exhumation along the Shuoqu River valley.

Data Dispersion

The standard deviations of the mean ages of the samples from the Cilincuo massif vary between 2 and 14 m.y. (23%) for the ZHe ages and between 0.8 and 5.3 m.y. (37%) for the AHe ages. Both the ZHe and AHe ages of the other two samples, S16X9 and S16X14, reproduced quite poorly, and S16X9 had the highest standard deviation for both sets of ages: 34.8 m.y. for the ZHe ages and 15.8 m.y. for the AHe ages. Age dispersion can reflect differences in the properties of the grains that affect the He diffusion kinetics, such as the grain size, the U and Th concentrations, and the grain composition. The grain size controls the volume of the diffusion domain such that within the same sample, larger grains have older ages than smaller grains (e.g., Reiners and Brandon, 2006). Differences in the U and Th concentrations result in variable amounts of accumulated radiation damage, which can cause the He diffusivity to decrease at moderate damage levels and to increase at high damage levels (Guenthner et al., 2013). Another source of age dispersion is the uncertainties in the alpha-ejection correction of the ages, which can be very large depending on the zonation of U and Th in the dated zircon (Hourigan et al., 2005). All these factors can affect independently and simultaneously the resulting ages, with amplified effects at low cooling rates, and it is often very difficult to discriminate the factor(s) that is(are) dominant.

Sample S16X14 is from a sandstone, and based on its Late Cretaceous ZHe ages, we speculate that the age dispersion of this sample was likely controlled by the combined effect of different grain composition and a long residence at low temperatures where radiation damage accumulated. Sample S16X19 is from a strongly weathered Triassic granite (site 4 in the DongWang area) with very low apatite quality, which largely explains the poor reproducibility of its AHe ages but not of its ZHe ages. Other samples from the same region and rock type or other data are not available to gain some insights into the scarce reproducibility of this sample, and therefore we do not further consider it for our interpretation.

The large number of grains dated from the Cilincuo massif allows a more detailed analysis of the possible causes of the age dispersion, which is, however, moderate for most samples. If the age dispersion in our samples is related to variable grain size, then we should observe mostly positive correlations between age and grain size. The equivalent spherical radius (Rs; e.g., Reiners and Brandon, 2006) of a grain is proportional to its diffusion domain, and in the dated grains, it varies between 30 and 80 μm, with a positive correlation with age in only two samples (Data Repository Fig. S2A). If our samples were mostly affected by radiation damage, then we should observe correlations between the ages and the effective U concentration (eU; Flowers et al., 2007), which is a proxy for radiative damage. These correlations should be positive at moderate eU values and negative at high eU values (Guenthner et al., 2013). All zircon grains had eU values between 130 and 1900 ppm, and they showed positive, negative, or no correlations with eU, independent from the magnitude of the eU content (Data Repository Fig. S2A). Thus, we cannot resolve whether either the effect of grain size or the radiation damage, or both, could explain the observed age dispersion (Data Repository Table ST3; Fig. S2). Some constraints on the possible extent of U and Th zonation of our samples can be derived from the comparison of the U and Th concentrations at the rim and at the center of the zircons from sample S16X10 from the Cilincuo massif that we dated with the U-Pb method. In this sample, the ratio between U or Th at the rim and U or Th at the center ranges from 0.5 to 5 for >90% of the measured zircons. For a fivefold increase of the U-Th concentrations at the rim, the alpha-ejection correction factor Ft calculated using the program of Gautheron and Tassan-Got (2010), and based on grain geometries similar to those of the zircons from sample S16X10 (in particular, S16X10-a, S16X10-b, and S16X10-c), would be 0.01–0.1 times smaller than the alpha-ejection correction factor calculated assuming no zonation. Thus, although we did not test the zonation for the other samples of the Cilincuo granite, we infer that zoning is not likely the main cause for the age dispersion of the Cilincuo samples. A possible explanation is, instead, that the large spread from 82 to 41 Ma of the ZHe ages from the Cilincuo massif may represent a long residence after intrusion of the Cilincuo granites at temperatures close to those of the ZHe PRZ. If the Cilincuo granites were intruded at shallow crustal levels, only a few kilometers deep, and if after intrusion they were affected by very slow exhumation, then they could have cooled rapidly down to temperatures close to those of the country rocks into which they intruded. In this scenario, the most retentive zircons would give ZHe ages close to the age of crystallization, whereas the least retentive zircons would record the later phase of slow cooling. In order to test whether this scenario is plausible, we performed thermal modeling using QTQt and we included in our models all the ZHe values that were not older than the crystallization age (thus, only two ages were discarded). Modeling results are discussed in the following section.

Thermal Modeling

The distribution of cooling ages within a landscape is largely controlled by variable exhumation rate coupled with a geothermal gradient and by differential exhumation represented by samples distributed across elevation (Wagner et al., 1979). An age-elevation correlation can be used to calculate apparent exhumation rates based on simple linear regression, or ages can be modeled for specific cooling histories using kinetic models for daughter product retention. The code used in this paper, QTQt, assumes that samples located closely together in elevation can be related through a constant temperature gradient and calculates ages using diffusion kinetics for the apatite and zircon systems. The dated samples from the Cilincuo pluton are located at elevations between 3300 and 3900 m (Figs. 7A and 7B), and all the ZHe ages vary between 80 and 40 Ma. At site 3, the linear regression through the mean ZHe ages of six samples distributed over 600 m of elevation offset has a low probability (R2 = 0.4) and a slope of 43 m/m.y. Thus, if we assume that, despite the low probability of the fit, the ages vary linearly with the elevation, then they would imply a very low exhumation rate. At site 1, the linear fit between the mean ZHe ages and the elevation has a negative slope. At site 2, the slope is positive with a moderate slope and with a poor probability. These observations suggest that if there is a correlation between cooling ages and elevation, it indicates a low exhumation rate. The AHe ages at elevations of ∼3900 m range between 40 and 30 Ma and at lower elevations center around 10 Ma. Unfortunately, like the ZHe samples, the AHe samples also show poor age-elevation correlations, even when only the ages around 10 Ma are considered. However, in contrast to the ZHe ages, the AHe ages below 3900 m elevation are distributed over a relatively narrow range that we interpret as indicating a phase of rapid cooling. Thus, we infer that the three sampled sites have similar cooling histories with a long residence at temperatures within the PRZ of the ZHe system during the Late Cretaceous–Paleogene followed by cooling below 100 °C starting at 40 Ma, with a rapid cooling phase during the Miocene or from the Miocene to the present. Some of the spread and misfit could be due to different retention characteristics of individual grains, e.g., radiation damage and/or grain size. We can account for this by modeling the data simultaneously using QTQt, as QTQt treats each grain differently using each grain’s specific information (size, isotopic concentrations, radiation damage). In Figure 7, we present our preferred modeling results, which were obtained by imposing starting temperature-time conditions (350 ± 100 °C; 75 ± 10 Ma) combined with a geothermal gradient (110 ± 40 °C) meant to simulate initial conditions possibly similar to those at the time of the intrusion. We purposely used a wide initial T-t range, and, except for site 3, we did not add other temperature-time constraints to the model because, as noted by Vermeesch and Tian (2014), QTQt requires correct assumptions on the model and is susceptible to overinterpretation, since it always produces results no matter what the initial constraints are. All three sites showed rapid cooling to temperatures below 200 °C during the Late Cretaceous following the emplacement of the granite, which we interpret as postintrusion cooling. For site 1, this episode was followed by a slow cooling and a prolonged stay at temperatures close to 100 °C. Thus, different from our initial inference, for the time after the rapid postintrusion cooling, this model predicts a prolonged stay at temperatures lower than, rather than within, the ZHe PRZ. This prediction probably reflects the facts that, at this site, all but one ZHe grain ages are older than 60 Ma, and the best model results fit the majority of data as shown in Figure 7B for site 1. In fact, in this figure, the predicted mean cooling path passes close to the oldest grain ages and does not fit well with the 45.5 Ma ZHe grain age of the highest sample. At site 2, the model indicates a long residence through the Paleogene to temperatures even lower than 100 °C. This is partly dictated by the ∼40 Ma AHe age of the highest sample at this site. However, we argue that given that eight out of 10 ZHe grain ages at site 2 lie in the range between 60 and 40 Ma, it is likely that cooling during the early-middle Paleogene occurred at temperatures between 200 °C and 100 °C rather than below 100 °C. Irrespective of the lowest temperatures possibly reached by our samples during the Paleogene, models indicate a rapid cooling shortly after intrusion followed by a long Paleogene stay at low temperatures that probably correspond to depths not greater than ∼6 km. This early phase was followed for all three sites by a rapid cooling episode during the middle Miocene or from the middle Miocene to the present. Unfortunately, our data and models can only poorly resolve whether the second fast cooling event occurred stepwise with a pulse around 10 Ma followed by slow cooling or continuously till the present with faster rates after 20 Ma. This is true especially for site 3, where several AHe grain ages are younger than 10 Ma. For this site, we imposed temperatures higher than 100 °C until 40 Ma to obtain a stepwise cooling as shown in Figure 7C; without this constraint, the model would predict continuous cooling, rather than a stepwise cooling path from the middle Miocene to the present. Thus, we speculate that if by the end of the Paleogene, the geothermal gradient in the Cilincuo pluton had equilibrated to close to a normal gradient, then it is more likely that the last cooling event occurred stepwise rather than continuously.

It is also important to note that for the thermal models presented in Figure 7, the initial geothermal gradient was allowed to vary only during the final cooling phase, resulting in an exaggerated offset among the samples throughout most of the cooling history. In order to know how much this initial condition influenced our results, we tested different parameters in additional models presented in the Data Repository Section 5. In particular, we performed thermal simulations where the temperature offset among the samples was allowed to vary at any time during the cooling history. We also tested models where the initial gradient was set at 30 ± 10 °C/km to simulate the fact that if the intrusion occurred at shallow depth, rapid cooling after the intrusion might have allowed the geothermal gradient to equilibrate very rapidly with the gradients of the host rocks.

Local Incision History and its Implications for the Cenozoic Exhumation History along the Yangtze River

From the thermochronometric and the sedimentologic data, we obtained two independent proxies for Cenozoic fluvial incision in the Shuoqu valley. Detrital sediments that partly fill a paleosurface 40 km south of the Cilincuo massif show variable incision of 0.5–1.2 km. Thermochronometric ages of the Cilincuo massif indicate a Miocene cooling event. If the exhumation of the Cilincuo massif was driven mostly by fluvial incision, then it could be related to the incision by the Shuoqu River, as constrained by the red beds. Thus, we could speculate that the exhumation of the Cilincuo massif constrains a maximum value for the fluvial incision, whereas the incision of the red beds represents a minimum estimate of the regional fluvial incision. However, there is no solid constraint on the age of these Cenozoic deposits, which prevents any quantitative comparison between these two close sites. Further field studies are necessary to unravel the deposition history of the red beds.

Based on AHe ages ranging from 9 to 7.2 Ma and Eocene–Oligocene ZHe ages, Ouimet et al. (2010) proposed that rapid incision by the upper Yangtze River started between 15 and 10 Ma and continued until present (Figs. 2 and 3). Based on a direct interpretation of their age-elevation relationship, Ouimet et al. (2010) argued for a constant incision rate of 0.3–0.4 km/m.y. from 5 Ma to the present. Thus, the Shuoqu River and the Yangtze main trunk seem to share a common exhumation history. If incision by both rivers was synchronous, it is possible that they were part of the same drainage network at the time of incision, i.e., by 15 Ma. However, this could also be due to a regional exhumation signal that affected the entire SE Tibetan Plateau synchronously, which implies that the Shuoqu and the Yangtze main trunk were simply incising synchronously without necessarily being connected to each other. We discuss the hypothesis of a synchronous exhumation in the next section.

Based on AHe ages ranging from 112 to 50 Ma, Clark et al. (2005) suggested that slow erosion occurred within the Yangtze basin downstream of the first bend (Fig. 3), where terranes are drained by tributaries of the middle Yangtze. These samples were collected at high elevations (>3860 m; Fig. 3A), which correspond to a relict low-relief surface (Clark et al., 2004) that was later dissected by river incision. Tian et al. (2014) sampled the same pluton as Clark et al. (2005) and obtained Cretaceous AHe ages for samples above 3700 m. They also collected two additional samples at lower elevation (3440–3150 m) in an adjacent valley. These two samples yielded AHe ages of 23 Ma and 15 Ma. It is likely that these younger ages record the onset of exhumation caused by river incision and that this event did not affect the higher-elevation samples of Clark et al. (2005), located in an area probably preserved from incision, as proposed by Clark et al. (2005). AHe data in that area present some similarities with our data (Fig. 3): Our low-elevation samples, located at elevations lower than 3900 m, record rapid Miocene cooling, initiated around 18–15 Ma, whereas higher-elevation samples reflect slow cooling.

If the upper Yangtze and the Shuoqu River were in fact connected by the middle Miocene, then the apparent difference in cooling histories during the last 5 m.y. recorded by our data and by those of Ouimet et al. (2010) could simply derive from the fact that the data from the two studies sampled different elevations. Most of our samples are from elevation >3300 m, whereas those of Ouimet et al. (2010) came from elevations ≥2300 and ≤ 200 m. Thus, we speculate that in fact the two data sets encompass different temporal and spatial intervals of the same history of incision in the upper Yangtze drainage: Incision started at ca. 15 Ma, affected the landscape that is now preserved along a tributary at elevations between 4000 and 3300 m prior to ca. 12 Ma, and it further entrenched along the main trunk from 8 Ma onwards.

Further south, McPhillips et al. (2016) proposed that incision of 1 km occurred near the Yangtze first bend in the Jianchuan area (Fig. 3) between 18 and 9 Ma, using burial ages of sediments collected in abandoned caves preserved above tributaries of the Yangtze River. These data are located close to the two bedrock AHe ages of Yang et al. (2016) of 38 Ma and 29 Ma and to the altitudinal profile of Shen et al. (2016), which yielded AHe ages that range from 31 to 18 Ma and have a positive correlation with elevation, interpreted as a cooling episode related to fluvial incision. Yang et al. (2016) also reported a single sample yielding a mean AHe age of 31 Ma (sample Yan08; Fig. 3) a few kilometers north. Last, Cao et al. (2019) collected two samples in the same area that yielded mean AHe ages of ca. 38 and 30 Ma. By modeling these ages along with fission-track data from an altitudinal profile, Cao et al. (2019) showed that a rapid cooling episode occurred from 28 to 20 Ma, leading to a total exhumation of 2.3–3.2 km. According to Cao et al. (2019), there has been less than 0.2 km of exhumation in the Jianchuan Basin since 20 Ma. If during the Miocene, incision was less than ∼2 km (as proposed by McPhillips et al., 2016), then it would have limited effects on the AHe system, which in turn would not record this late incision phase. This could explain why the thermochronometry data (Yang et al., 2016; Shen et al., 2016; Cao et al., 2019) predate the incision history recorded by the burial dating conducted by McPhillips et al. (2016). Thus, it is possible that at least two exhumation episodes occurred near the first bend between the Eocene and the Miocene and that uplift (or at least exhumation) occurred through distinct phases.

Implications on the Exhumation History of SE Tibet

Although there has been an increase in thermochronometric studies in SE Tibet, the data set is still scarce with respect to the region’s size, and the synchronicity of exhumation and the relative importance of tectonic versus climatic forcing in controlling the Cenozoic exhumation are debated. Part of the thermochronometric data set is clearly related to active tectonics. For example, the exhumation along the Litang fault system accelerated during the late Miocene–early Pliocene because of normal motion on segments of the fault (Figs. 2 and 3; Zhang et al., 2015). Clark et al. (2005) proposed that rapid river incision initiated during the middle Miocene simultaneously throughout SE Tibet, as a result on the uplift of that region of the plateau. Clark et al. (2005) suggested that resulting higher elevations were related to the increase in SE Asia monsoon intensity, which has been attested by weathering records in the Indus and Bengal Fans (Clift et al., 2008). Nie et al. (2018) also supported the hypothesis of a link between monsoon intensification and rapid river incision recorded by AHe data along the Mekong River; however, this was based on the comparison of one AHe data point in the upper reaches of the Mekong River with six samples from an altitudinal profile located ∼800 km further south along the lower Mekong River (Fig. 3A). In Nie et al.’s scenario, increased monsoonal precipitation during the middle Miocene triggered a simultaneous response along the entire course of the Mekong River, leading to a lowering of the river longitudinal profile that the authors modeled using a stream power law. Because Asian monsoon presently affects most of the SE Tibetan Plateau, enhanced monsoon during the middle Miocene should have operated similarly at the regional scale. In that case, one could expect to observe a similar timing of exhumation in the Salween, Mekong, and Yangtze Rivers. However, available data, including Pliocene AHe ages in the Salween area (Yang et al., 2016) and an Oligocene exhumation phase near the Yangtze in the Jianchuan area (Shen et al., 2016; Yang et al., 2016; Cao et al., 2019), rather point to a heterogeneous pattern of exhumation that is hard to reconcile with a climate driver. In the case of the Jianchuan Basin, Cao et al. (2019) suggested that late Oligocene–early Miocene exhumation was due to approximately N-S–oriented thrust faults east of the basin (Fig. 2). In the Mekong area, Liu-Zeng et al. (2018) and Yang et al. (2016) combined AHe, ZHe, and apatite fission-track dating that indicate two exhumation phases ca. 60–40 Ma and ca. 20 Ma to present. The picture emerging from the increasing number of thermochronometric studies that improve the spatial coverage of SE Tibet points to an apparent heterogeneity in space and time of the Cenozoic exhumation of the southeastern Tibetan Plateau margin.

Cenozoic sediments in the Xiangcheng area consist of sandstone and unsorted conglomerates that are locally sourced and that unconformably rest on a paleosurface. Postdeposition incision of these sediments varies between 0.5 and 1.2 km. Thermal modeling of bedrock apatite and zircon (U-Th)/He thermochronometry near Xiangcheng city, SE Tibetan Plateau, shows that this region experienced a pulse of exhumation that initiated around 18–15 Ma and ended ca. 12 Ma. The Xiangcheng area is affected by a local sinistral strike-slip fault that does not influence the local exhumation pattern. The middle Miocene exhumation phase revealed by our results was due to fluvial incision by the Shuoqu River. We suggest that the upper Yangtze and its tributary (Shuoqu) were already connected by 15 Ma. At the regional scale, a comparison of our results and published thermochronometric studies suggests that the Cenozoic exhumation in the southeastern margin of the Tibetan Plateau is heterogeneous in space and/or time.

This study was funded by ETH postdoctoral fellowship FEL-11 15–2 (to Gourbet), the National Natural Science Foundation of China (41602210, 41472182), and the Fundamental Research Funds for the Central Universities (2018XZZX001–03). We thank Lu Gang (ETH), Wang Yanyan (ETH), Wang Xinsong (McGill University), and Guangchun Fei (Chengdu University of Technology) for their help with the Chinese literature, and Marie-Luce Chevalier, Anne Replumaz, and Sean Gallen for helpful discussions.

1GSA Data Repository Item 2019409, field observations and samples locations, isotopic data for zircon U-Pb dating, an assessment of error and data dispersion of (U-Th)/He ages, and additional thermal models discussed in the main text, is available at http://www.geosociety.org/datarepository/2019, or on request from editing@geosociety.org.
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.