To reconstruct the complete exhumation history of the SE Tibetan Plateau, we compiled 1,202 low-temperature thermochronological data from previous studies throughout the study area, including AHe (n = 260), ZHe (n = 182), AFT (n = 633), and ZFT (n = 127) (Figures 1B, 2). All of these ages are considerably younger than the corresponding stratigraphic or crystallization ages, which is a reflection of the post-depositional or monotonous cooling histories of the basement rocks. It should be noted that we only consider the Cenozoic cooling and exhumation of the SE Tibetan Plateau. Ages older than 64 Ma are not employed for the qualitative investigation of cooling age patterns and inverse modeling as they provide little information (e.g., Herman et al., 2013; Fox et al., 2016; Jiao et al., 2017; Willett et al., 2020). Furthermore, we split the area into seven major terranes/units based on geological formations (Figure 1A; Burchfiel and Chen, 2012) and sample distribution, which include Songpan–Ganze–Yidun, Longmen Shan–Sichuan Basin, Qiangtang–Lhasa, transitional zone, Kang-Dian, Western Yunnan, and Yangtze (Figure 1B). Of these, we have not collected low-temperature thermochronological data for the Yangtze, so we focus our analysis on the other six tectonic units.

KDE is a non-parametric method for calculating the probability density function of random variables. KDE is a continuous and smooth substitute for the standard histogram. KDE is a function that uses intrinsic adaptive bandwidth to account for analytical uncertainty and overlays a Gaussian “bell curve” on top of each measurement. These adaptive functions eliminate the uncertainties generated by analytical errors and age abundance by adjusting the bandwidth dependent on the local density. Vermeesch (2012) provided a detailed explanation about KDE. Therefore, it will not be repeated in this paper.

For the quantitative assessment of the spatiotemporal exhumation pattern of the SE Tibetan Plateau, we employed a Gaussian linear inversion approach based on compiled AHe, ZHe, AFT, and ZFT ages. This inversion approach was established by Fox et al. (2014). We adopted the updated version described in Herman and Brandon (2015). This approach exploits the exhumation rates from age-elevation sections, as well as multiple thermochronological systems with different closure temperatures (Fox et al., 2014). According to this approach, depth (Z _c ), from the closure isotherm/depth to the surface, is defined as the integral of exhumation rates (ė) over a finite number of time intervals, from the cooling age (τ) to the present day. The equation is represented as follows: Z c = ∫ 0 τ e ˙ d t . (1)

For providing a numerical solution for Eq. 1, we need to discretize the time integral before parameterizing the exhumation rate as a piecewise constant function. The subsequent is as follows: ∑ i = 1 N − 1 e ˙ j ∆ t + e ˙ N R = Z c , (2)

where e ˙ j represents the average exhumation rate within uniform time intervals (∆t). N is the number of intervals, R is the remainder of the division of τ by N-1, and R = τ-(N-1)∆t.

However, the model assumes that the rocks experience monotonic cooling histories and do not account for the reburial process (i.e., negative exhumation rates). Hence, Eq. 1 could be converted to a logarithmic space to prevent this situation: ζ = ln ∫ 0 τ ⁡ exp ε d t , (3)

where ζ = ln Z c and ε = ln e ˙ .

The parameter e ˙ or ε in the equations is an unknown quantity, and solving the linear inversion is required. Therefore, we provide a priori exhumation rate ( e ˙ p r i ) and obtain an estimate of the variance about this mean (σ²). Hence, the a priori model covariance matrix is built using this variance: C M i , j = σ 2 ⁡ exp ⁡ ⁡ − d L 2 . (4)

The matrix C M i , j is constructed for each time interval containing the horizontal distance (d) between the ith and jth samples, and L represents a specified correlation length.

Finally, a non-linear least square method is used to solve the exhumation rate by employing multiple iterations: ε m + 1 = ε m − μ C M G t C D − 1 ζ m − ζ o b s + ε m − ε p r i , (5)

where m is the iteration number, ε p r i is the logarithm of the initial value as the a priori exhumation rate ( e ˙ p r i ), G is the product of the exponential of ε and time interval divided by the closure depth (Herman and Brandon, 2015), C_D is the data covariance matrix, μ is an arbitrarily chosen parameter which controls how rapidly the model parameters change, and ζ m and ζ o b s are the logarithms of the modeled and observed closure depths, respectively.

As the exhumation history obtained from the inversion may depend on the chosen priori exhumation rate, we used various values of the parameters, such as its variance, the initial geothermal gradient, the spatial correlation length, and the time intervals, to test the robustness of our results (Supplementary Table S1). A total of five inverse models are depicted in Supplementary Table S2, which is provided as a supplementary material. The M₁ model, which provides the optimal solution out of these models, is the one on the basis of which the results have been interpreted. The detailed interpretation is described in Section 5.1. In the M₁ model, we used a prior exhumation rate of 0.1 ± 0.1 km/Myr, a correlation length scale of 20 km (Yang et al., 2016), and a time interval length of 4 Myrs, from 64 Ma to the present. The tested prior exhumation rate is the end-member value for the regions that exhibit higher exhumation rates like the Jianchuan Basin (0.06 km/Myr before Early Cenozoic; Cao et al., 2020). Both a spatial correlation length and time interval act as trade-off parameters between solution resolution in space and the need to smooth the noise (Yang et al., 2016). A larger correlation length will always increase the numerical values of the resolution and decrease the variance, whereas the time interval number will always decrease the resolution (Willett et al., 2020). Therefore, it is necessary to choose the aforementioned suitable values for simulation. In addition, we also set a fixed model thickness of 30 km with a temperature value of 10°C at the surface (Liu et al., 2021; Shen et al., 2021), which produces an initial geothermal gradient of 30°C/km (Yang et al., 2016). The thermal diffusivity is 25 km²/Myr (Braun et al., 2012), and heat production is not considered; as reported by Batt and Braun (1997), the curvature of the geotherm caused by heat-producing elements can be neglected while interpreting low-temperature thermochronology.

The 1,202 available low-temperature thermochronological data were statistically analyzed using KDE with various terranes/units in order to clarify the full exhumation signals of the SE Tibetan Plateau. This analysis ultimately produced a multi-peak distribution of age patterns (Figure 3). Of these, a total of 307 data pertained to the Songpan–Ganzi–Yidun Terrane, revealing five peaks (Figure 3A): 61 Ma (ZHe), 32 Ma (AFT/ZHe), 13 Ma (AFT/ZHe), 8–6 Ma (ZFT/AHe), and 4 Ma (AFT). Furthermore, 478 data pertained to the Longmen Shan–Sichuan Basin, with four population peaks since the Oligocene (Figure 3B), including 27–23 Ma (AHe/ZHe), 17–13 Ma (ZFT/AHe), 10–8 Ma (AHe/ZHe), and 3 Ma (AFT/AHe). With respect to the Qiangtang–Lhasa Terrane, AHe data revealed two Miocene age peaks at 18 Ma and 8 Ma, while the AFT ages revealed one population cluster at 4 Ma (Figure 3C). In contrast, the sample size of the transitional zone is even smaller, with only 13 data (AFT = 10 and AHe = 3); thus, only the AFT data were statistically analyzed, and based on previous studies, it was concluded that two peaks existed in this area, at ∼15 Ma and ∼4 Ma (Clark et al., 2005; Lai et al., 2007; Zou et al., 2014; Figure 3D). Furthermore, 46 data pertaining to the Kang-Dian Terrane (AFT = 22, AHe = 15, and ZHe = 9) were available; on analyzing the data, four peaks were revealed (Figure 3E): 36–35 Ma (AFT/ZHe), 27–24 Ma (AFT/ZHe), 18–15 Ma (AFT/AHe/ZHe), and 11–6 Ma (AFT/AHe/ZHe). As compared to other terranes/units, the western Yunnan Terrane generated six age peaks. The oldest peak was recorded by the ZFT ages at ∼58 Ma, which was followed by the AFT, AHe, and ZHe age clusters at 38 Ma. In addition, four other young peaks were reported, which included 30–24 Ma (AFT/AHe), 13 Ma (AHe/ZHe), 6 Ma (AFT), and 4 Ma (AHe) (Figure 3F).

In summary, the results revealed synchronous exhumation events in the different terranes/units of the SE Tibetan Plateau despite the existence of variable tectonic environments. We, therefore, suggest that the SE Tibetan Plateau might have experienced at least six rapid exhumation events since the Cenozoic, including ∼61–58 Ma, 38–35 Ma, 32–23 Ma, 18–13 Ma, 11–6 Ma, and 4–3 Ma.

The results presented in Figures 4–9 represent one solution (M₁) amongst a series of acceptable solutions utilizing various parameters or a priori models. However, previous studies have suggested that some factors, such as the geothermal gradient, prior exhumation rate, time interval, and data search radius (correlation length), could affect the modeled results (Fox et al., 2014; Ballato et al., 2015; Yang et al., 2016; Siravo et al., 2019; Stalder et al., 2020). To verify the reliability of our inversion, we performed a sensitivity study to explore the aforementioned parameters (Supplementary Tables S1, S2; Supplementary Figures S1–S4).

Figures 10, 11; Supplementary Materials indicate that different parameters may affect the inversion results. First, an increase in a geothermal gradient from 25°C/km to 30°C/km may lead to a lower exhumation rate because heat is more efficiently advected to the surface, and thus, the closure depths are shallower (Figures 10A, D). Second, the variations in the a priori exhumation rate can potentially be relevant. In case of inadequate data resolution (fewer data; Figures 11A, B), these data are not sufficient to restrict the exhumation rate well within the function interval (sampling interval) during the calculation. Consequently, the model solution is particularly sensitive to the a priori model (Figures 10A, B). In this case, our model values do not deviate from the a priori parameter. In other words, at the initial stage of calculation, when the prior exhumation rate was large, the calculated result was larger than the small one. Alternatively, as the resolution of the data increases, there may be an artifact of a sharp increase or decrease in the exhumation rate (Fox et al., 2014; Willett et al., 2020). In addition, the time interval and the correlation length of the samples also affected the modeled results. Smaller time intervals and shorter sample correlation lengths reduce the amount of data available to calculate exhumation rates in a given area, increasing the value of data variance and reducing resolution. On the contrary, when both of them increase, restriction of exhumation rates in the calculated area was facilitated (Fox et al., 2016; Stalder et al., 2020; Figures 10A, C).

However, the aforementioned results only have a strong influence on areas with lower resolution (Figures 10A–E, Figures 11A–E) but are weaker for the higher-resolution areas (Figures 10F–J, Figures 11F–J). Furthermore, previous studies have shown that a priori values have little impact on the final simulation results (Fox et al., 2014; Willett et al., 2020). To ensure reliability of the results, we only analyzed the areas having more than 0.25 resolution (Herman et al., 2013).

Similarly, we have illustrated the frequency of the difference between predicted and measured ages (misfit) for the different models through a histogram (Figure 12). Notably, the results pertaining to different parameter simulations did not vary obviously (Figure 12), and the differences mainly ranged between −5 Ma and 5 Ma. We statistically analyzed the misfit values of the various models independently, to examine the impact of each parameter on the results pertaining to age prediction, and the outcomes; the results are illustrated in Supplementary Table S2. M₁, which considers the data dispersion and practical needs, was ultimately selected as the optimal model for analysis, even though the findings pertaining to other parameters were slightly different.

Our findings suggest that the rapid Cenozoic cooling episode initially occurred during the Paleocene (Figure 3). This event was recorded not only in the SE Tibetan Plateau but also in the central, northern, and eastern parts of the plateau (Hetzel et al., 2011; Rohrmann et al., 2012; Dai et al., 2013; Haider et al., 2013; Jian et al., 2018). Such a synchronous and widespread exhumation suggested a regional event at this stage: crustal shortening due to large-scale tectonic activity with the beginning of the collision between the Indian and Eurasian continents (England and Searle, 1986; Murphy et al., 1997; Kapp et al., 2005; Kapp et al., 2007). Although the onset of the India–Eurasia collision remains controversial, the prevailing viewpoint suggests 65–55 Ma (Wang et al., 2014), which overlaps with the Paleocene rapid exhumation event. During this period, the compressive stress trending N–S was promptly transmitted to the northern margin of the plateau, activating the former and deep faults, leading to thrust faulting and folding, and characterization of the orogen, as illustrated in the low-temperature thermochronology data (e.g., He et al., 2017; Liu and Cai, 2022). These views support the immediate response model for growth of the Tibetan Plateau (e.g., Dayem et al., 2009) but contradict the models that support diachronous progressive growth from south to north (e.g., Tapponnier et al., 2001).

By the Early–Middle Eocene, our KDE results could not depict an age peak, but higher exhumation rates were observed in the local regions (i.e., Lincang and Xiangcheng). Furthermore, such a phase of exhumation was also reported at Weixi and Deqin (Figure 13A), with a rapid exhumation rate of ∼0.1–0.3 km/Myr, as revealed by AFT and AHe data (Liu-Zeng et al., 2018b). It is worth noting that Ou et al. (2020) recalculated a rate of 0.04–0.06 km/Myr, at this stage, using 3D thermo-kinematic modeling near Deqin, and their findings are consistent with our findings obtained from inversions (Figure 5; Supplementary Materials). Nevertheless, the report of the Early–Middle Eocene phase of exhumation remains limited and suggests that it is probably a local signal (Zhang et al., 2022).

During the Late Eocene–Early Oligocene, the linear inversion results revealed that the Xichang Basin and Lincang generated enhanced exhumation rates (Figure 6), which were confirmed by previous results and our KDE results (Figures 3, 13A). We concluded that a tectonic shortening mechanism may explain why the aforementioned two regions experienced rapid uplift during this phase (Deng et al., 2018; Liu et al., 2021). At this stage, the Indian plate continued to move northward, and under the influence of oblique compressional collision, strong tectonic deformation began in the SE Tibetan Plateau, subsequently leading to development of numerous thrust-nappe tectonic belts and strike-slip shear zones, accompanied by a series of contemporaneous igneous rocks (Wang and Burchfiel, 1997; Tapponnier et al., 2001). These phenomena have been observed in Lincang and Xichang Basins. For example, the low-temperature thermochronological data identified that the Xichang Basin underwent protracted cooling at 40–20 Ma with a rate of 0.1–0.3 km/Myr (Deng et al., 2018), and the findings are consistent with our results (Figures 6D–F). Meanwhile, they also recovered the stratigraphy from the Late Triassic in and around the Xichang Basin by balanced cross-sections and hypothesized some shortening of ∼20–30 km occurred in the east–west direction. Furthermore, there are only two unconformities among the upper Triassic and pre-Late Triassic sequences (T₃/P), Neogene, and its underlying strata (N₂/E₁) in the Xichang Basin (Deng et al., 2018). Considering the aforementioned points, we speculated that the shortening should be explained in terms of Cenozoic shortening, rather than an average shortening since the Late Triassic. Meanwhile, the intermontane Lanping–Simao Basin was folded, closed, and converted into a mountain, with simultaneous molasse sediment signals in the front of the mountains between the Late Eocene and Oligocene (He et al., 1996; Wang and Burchfiel, 1997). In addition, the abundant Eocene–Oligocene thrust-nappe tectonic belts were developed in the interior and margin of the Lincang granite batholith, which was composed of several thrust or oblique faults with right strike-slip characteristics (Li, 1994; Wu, 1994; Liu et al., 1998). Amidst such thrust-nappe activities, the Lincang granite batholith was uplifted and exhumed at a large scale (Shi et al., 2006). Similarly, upper Eocene to Oligocene detrital strata overlying the older deformed strata underwent regional unconformity in the Simao and Baoshan/Lincang blocks of the SE Tibetan Plateau, indicating the beginning of shortening deformation (Wang and Burchfiel, 1997).

The Oligocene cooling phase is coeval with those reported around the SE Tibetan Plateau, suggestive of enhanced exhumation on a regional scale (Figure 3). At the northeastern terminus of the Nantinghe fault, this Oligocene cooling was determined by thermal modeling based on the AHe age (Wang et al., 2019; Sample LC03-1), which may result due to tectonic transtension/transpression toward the northeastern termination. Cao et al. (2020) revealed rapid exhumation of 2.3–2.5 km at ∼28–20 Ma (Figure 13A) in the north of the Jianchuan Basin. They inferred the current relief of the SE Tibetan Plateau can be triggered by the Yulong thrust belt, while the possibility of crustal thickening and increase in elevation before the Oligocene or Early Cenozoic collision of India–Asia cannot be ignored. In contrast, Shen et al. (2016) recognized a phase of exhumation in the Jianchuan Basin between 30 and 20 Ma with a rate range of ∼0.10–0.18 km/Myr (Figure 13A), resulting from the transpressional deformation in the SE Tibetan Plateau.

The Early Oligocene rapid exhumation and uplift of the Lincang granite batholith overlay with the dextral/sinistral strike-slip movement, such as Gaoligong, Chongshan, and Ailao Shan–Red River shear zones, even though the initial timing of their activities is debatable (Leloup et al., 1993; Leloup et al., 2001; Gilley et al., 2003; Akciz et al., 2008). Since the Oligocene, the rocks located in these shear zones underwent shortening, extruding toward the southeast, and differential clockwise rotation, consistent with the current crustal motion derived by GPS (Pan and Shen, 2017). Such kinematic properties indicate that they are significant to absorb hundreds of kilometers of displacement of the SE Tibetan Plateau (Wang and Burchfiel, 1997; Yin and Harrison, 2000; Gilley et al., 2003; Akciz et al., 2008).

We observe that since the Middle Miocene, the entire SE Tibetan Plateau exhibited significantly accelerated exhumation, particularly in the Gaoligong shear zone, the middle segment of the Longmen Shan fault, and the Xianshui River fault (Figure 8). Ge et al. (2020) used thermal history inversion to derive a rapid exhumation rate of 0.4–0.6 km/Myr at this stage, and our simulation results align with the finding (Figure 8). They suggested that such an enhanced exhumation phase should be attributed to the tectonics, combined with previous ⁴⁰Ar/³⁹Ar and U-Pb dating (Akciz et al., 2008; Lin et al., 2009). However, we concluded that this phase of rapid cooling was caused by intensified monsoon precipitation (Nie et al., 2018; Liu et al., 2021). At 17–14 Ma, the globe was in a period of Middle Miocene Climate Optimum, and increased CO₂ levels and warm climate were the main factors contributing to enhanced precipitation from the Asian summer monsoon during this period (Licht et al., 2014; Ji et al., 2018). Enhanced precipitation will potentially increase the power of river incision and trigger fault activity, which itself can intensify incision by increasing relief (Eroğlu et al., 2013; Wang et al., 2019). This hypothesis has been clarified by Ren X. et al. (2020). For the Xianshui River fault, previous work has dated U–Pb, mica K–Ar, and AFT to indicate that it initiated to undergo sinistral strike-slip movement at ∼15–13 Ma, suggestive of the onset of compression toward the western Sichuan and South China Block (Roger et al., 1995; Wang et al., 2009).

In the Late Miocene, the exhumation rates increased throughout the SE Tibetan Plateau, as compared to the previous phase, particularly in areas lying on the tectonic boundary areas (Figure 9). In addition, the analysis in Figure 3 shows that the entire SE Tibetan Plateau is in a phase of rapid exhumation, which was more widespread than that in the Middle Miocene. Similarly, the widespread rapid exhumation events coincide with the record of intense summer monsoon precipitation in Asia (Ren X. et al., 2020). The increased precipitation events of the Late Miocene may be attributed to an enhanced monsoon driven by tectonic uplift along the NE Tibetan Plateau compared to the enhanced monsoon brought about by global warming in the Miocene (Prell and Kutzbach, 1992; An et al., 2001). Topographic uplift helped in overcoming the effects of a cooler climate and decreased moisture level in the air (Ji et al., 2018) by intensifying the sea–land pressure gradient (An et al., 2001). Therefore, it was hypothesized that the uplift-driven monsoon intensification may lead to the rapid incision of the entire SE Tibetan Plateau. Furthermore, we collected paleoclimatic data for the Late Miocene period in Yunnan, when higher temperature and precipitation occurred as compared to the present era, with an annual mean precipitation of approximately 1,000–1,600 mm, which is indicative of a strong precipitation effect in Yunnan during the Late Miocene (Figures 13B, C). However, the tectonic possibilities cannot be ruled out. For example, in the Gaoligong–Chongshan shear zone, Wang et al. (2019) used apatite and zircon (U-Th)/He to derive two phases of rapid exhumation history in the area since the Late Cenozoic, one of which was prior to 11 Ma (15–11 Ma) and the other at 8–6 Ma (Figure 13A). The earlier one is linked with ductile transpressional shearing. However, the exhumation event at the Dayingjiang fault was caused by the sinistral transtension of a strike-slip fault (Wang et al., 2019). In the farther northern areas of the SE Tibetan Plateau, widespread rapid exhumation events were witnessed during the Late Miocene (Figure 13A), such as in the Longmen Shan (Kirby et al., 2002; Godard et al., 2009), Jiulong (Zhang et al., 2016), Min River (Tian et al., 2015), Min Shan (Tian et al., 2018), Xianshuihe fault (Zhang et al., 2017), Maoergai fault (Ansberque et al., 2018), Dadu River (Clark et al., 2005; Ouimet et al., 2010; Tian et al., 2015), the Daliang Mountains (Deng et al., 2014), and Anning River (Clark et al., 2005). Tian et al. (2013, 2015) suggested that time synchronization between the large-scale exhumation and rock and surface uplift in the Late Miocene was closely associated with crustal shortening rather than climatic effect, highlighting the role of crustal shortening in forming the enhanced erosion and high elevation of the plateau margin. In contrast, Clark et al. (2005) inferred that this synchronous event supported the tectonic (lower crustal flow)-climate (Asian summer monsoon) models.

During the Pliocene, despite the absence of significant fault activities in the Gaoligong–Chongshan shear zone, a rapid cooling event occurred in the region, which is thought to have been forced by climate-driven rapid incision of rivers or the injection of lower crustal flows (Wang et al., 2018; Ge et al., 2020). In contrast, Eroğlu et al. (2013) suggested that this rapid exhumation was triggered by crustal root delamination and the opening of the Andaman Sea. Additionally, in the Deqin–Zhongdian area, Xiao et al. (2015) and Zou et al. (2014) used AHe and AFT dating to derive a rapid cooling event since the Pliocene (5 Ma), and the calculated exhumation rates lay in the range of 0.3–0.4 km/Myr; our simulation results are consistent with the findings (Figure 9). However, there is an obvious difference in the explanation provided. Xiao et al. (2015) attributed the event to intense erosion due to strong precipitation from the Indian monsoon and the onset of the global glacial period. In contrast, Zou et al. (2014) revealed that it was caused by tectonic uplift, but no explanatory reason has been provided. In the Longmen Shan fault, the eastern margin of the Tibetan Plateau underwent a phase of eastward expansion due to strike-slip fault deformations and thrusting, and upper crustal deformation resulted in the formation of the Ganzi–Litang margin in the Miocene and continued to move toward the Longmen Shan fault around 5 Ma (Lai et al., 2007). The Xianshuihe fault, which intersects the Longmen Shan fault, played an important role in the tectonic evolution of the Tibetan Plateau. As discussed previously, the Xianshuihe fault began to undergo sinistral strike-slip movement at 15–13 Ma and was reactivated at 5 Ma, thereby cutting through the Longmen Shan fault and extending southward into the Xiaojiang fault zone (Roger et al., 1995; Wang et al., 2009).

Since 4 Ma, exhumation rates have experienced a general increase across the SE Tibetan Plateau (Figures 9D–F). An et al. (2008) suggested that the event was controlled by the tectonics, while other authors suggested that climate was the main driving force that controlled the river incision but did not exclude a tectonic role. Herman et al. (2013) used the Glide software application to invert the global cooling event by collecting nearly 18,000 low-temperature thermochronology data from different regions of the world. They found that erosion rates increased in all regions from 6 Ma, with a particularly noticeable increase from 2 Ma, and suggested that glaciation played a crucial role in this global increase in erosion rates. In the Qaidam Basin at the NE Tibetan Plateau, Ren X. et al. (2020) synthesized previous studies and demonstrated that a total of two intense precipitation events occurred at 4–2.7 Ma and 0.4 Ma. They illustrated that the former event was attributed to the further uplift of the Tibetan Plateau or closure of the Panama Seaway, whereas the latter may reflect a weak El Niño-like phenomenon related to high zonal sea surface temperature gradients. Therefore, we argue that climate forcing can play a key role in this cooling episode, along with minor tectonism.