Uplift History of the Eastern Pamir Inferred from Inversion of Thermochronometric Data and River Profile

The Pamir salient constitutes the northwestern indenter of the India-Eurasia collisional zone, forming an integral part of the Himalaya-Tibetan orogen [1, 2] (Figure 1). As one of the most tectonically active regions in Central Asia, the Pamir has accommodated a great amount of India-Eurasia convergence during the Cenozoic, leading to significant crustal deformations and complex deep geodynamic processes (e.g. [3, 4]). The Pamir may have experienced ~600–900 km crustal shortening during the Cenozoic [1, 2], whereas some recent observations revealed ~100 km of upper crustal shortening [5, 6] and ~50 km [5] to ~80–110 km [7] of northward translation since the Paleocene. Several large-scale domes consisting of crystalline basement crop out across the Pamir (e.g. [6, 8-11]; Figure 1). Metamorphic petrology and thermochronology studies have been conducted on the high-grade metamorphic rocks which document crustal thickening and peak prograde metamorphism followed by subsequent dome exhumation from a depth of ~30–40 km to the surface since Oligo-Miocene times (e.g. [2, 9, 12-14]). Current crustal deformation in the Pamir is characterized by thrusting and strike-slip faulting along its northern and western flanks and normal and strike-slip faulting in its interior and eastern margin [4, 15-19] (Figure 1). In contrast to the Tibetan Plateau with vast internal drainage basins of high elevation and relatively low internal relief, the Pamir is mostly externally drained, characterized by high relief and deeply dissected landscapes [20-23].

With its complex tectonic history, currently seismically active setting, and prominent topography, the Pamir salient is an exceptional site for studying Cenozoic orogenic process [24-26]. Over the past few decades, a wealth of geochronological, thermochronological, and structural data has been reported to reconstruct tectonic evolution of the Pamir (e.g. [6, 8-13, 20, 27-31]). In this study, we focus on the Eastern Pamir, where the crustal deformation over the late Cenozoic is dominantly controlled by the east-west extension along the Kongur Shan extensional system (KES), accompanied by exhumation of Kongur Shan and Muztagh Ata gneiss domes in its footwall [8, 9, 12, 20, 27, 28, 32, 33]. These structures make up the most prominent features, for example, a graben over 250 km long (i.e. the Muji-Tashkurgan Valley) and a population of peaks >7000 m a.s.l., within the eastern flank of the orogen. Several models have been proposed to elucidate the kinematic process that led to such tectonic and topographic anomalies within a convergent tectonic regime. The models of radial thrusting [16], synorogenic extension [32], and oroclinal bending [34] associate the interior extension in the Eastern Pamir with the motion of the Main Pamir Thrust (MPT) [16, 32] or the oroclinal bending of the entire orogen [34]. By comparison, some models proposed that the extension is due to northward propagation of the right-slip Karakoram fault [16, 35], while some researchers argued against such kinematic linkage between the strike-slip fault and extensional system [36-38]. Although a line of studies have provided estimates on regional tectonics in terms of timing, magnitude, and movement/exhumation rates (e.g. [8, 9, 14, 16, 20, 28-30, 33, 37, 39-42]), further research is required to establish a more detailed history of tectonic evolution in the eastern flank of the Pamir.

To better understand how current features within the Eastern Pamir developed, we report new thermochronometric ages along the Gez River transect across the Eastern Pamir. These new data, together with the published thermochronological ages (Online Supplementary Tables S1 and S2), were integrated with the longitudinal river profile to invert for regional uplift history along the Gez River from the KES to the MPT (boundary between the Pamir and the Tarim). Thus, our results provide a quantitative long-term estimate for the rock uplift rates across the Eastern Pamir.

The most prominent feature in the Eastern Pamir is the roughly north-south striking KES, which is ~250 km long and separates the wide and flat Muji-Tashkurgan Valley at ~3000–4000m in its hanging wall from the high mountains at >5000 m in its footwall. This extensional system consists of the NWW-SEE-striking dextral-normal Muji fault to the northwest, the west-dipping Kongur Shan normal fault and Tagharma normal fault, and the southernmost east-dipping Tashkurgan normal fault. Two gneiss domes, the Kongur Shan dome in the north and the Muztagh Ata dome in the south, are exposed in the footwall of the Kongur Shan normal fault defining the locally highest mountains (the 7649m Kongur Tagh [Kongur Shan] and the 7546m Muztagh Ata). The Eastern Pamir is mainly dissected by two rivers, the Tashkurgan River to the south and the Gez (Ghez) River to the north. Both rivers have their headwaters in the Muji-Tashkurgan Valley and then cut across the footwall of the Kongur Shan normal fault in their middle reaches. The riverbed reaches elevations of <2000 m as it enters the Tarim Basin. Previous studies suggest that the Gez River has been an antecedent river since the Oligocene [43, 44].

As the Gez River gorge has served as the main access to the Eastern Pamir, abundant field observations and thermochronological and petrological analyses were carried out along this transect in past decades (e.g. [8, 14, 20, 29, 31, 32]). Along the Gez River, a series of primary structures, including the Kongur Shan fault, Kongur Shan dome, Gez (Ghez) fault, Oytagh (Oytag/Wuyitake) fault, and the MPT from west to east can be observed (Figure 2). The Th-Pb dating on monazite grains from the metamorphic rocks of the Kongur Shan dome suggests peak metamorphism of ~8 kbar and 650–700°C at ~9 Ma [9]. Exhumation of the Kongur Shan dome is coeval with the east-west extension along the Kongur Shan fault. Yet, researchers have not reached a full consensus about the initiation of slip on the Kongur Shan fault due to limited records and different interpretations of the chronological data [8, 9, 20]. Fast cooling starting at 8–7 Ma recorded by K-feldspar ⁴⁰Ar/³⁹Ar sample from the Kingata Tagh massif is interpreted as the onset of the extension [9, 29]. Together with ~6–5 Ma cooling ages in the west flank of the Muztagh Ata dome (i.e. the footwall of southern Kongur Shan fault) [8, 27, 28, 30], it might imply a southward propagation of the initiation of east-west extension (e.g. [27]). In contrast, a dominant detrital zircon fission-track (ZFT) age group of ~6–4 Ma in fluvial sediments from rivers draining the gneiss domes provides evidence for a contemporaneous onset of the extension along the Kongur Shan fault at 6–5 Ma [8, 21]. A sharp increase in the proportion of detrital ZFT ages >10 Ma further to the north along Kingata Tagh massif demonstrates that the Kalagile fault forms the northern boundary of the rapid exhumation dome structure [8] (Figure 2). Given that the rapid cooling at 8–7 Ma [9] is located north of the Kalagile fault and similar cooling ages of 10–7 Ma occur at thrust zone within the Eastern Pamir [8, 30], this rapid cooling at 8–7 Ma may not necessarily date initiation of the Kongur Shan fault. Rather, it could be related to local exhumation associated with Late-Miocene thrusting [8].

Eighteen new bedrock samples were collected between the Kongur Shan normal fault and the MPT along the Gez River. These samples are distributed over a horizontal distance of ~60 km (Figure 2). Among these samples, nine (GZ1601–GZ1609) were from a vertical transect with an elevation difference of ~800 m, and the rest were collected at elevations close to the thalweg of the Gez River.

These samples were processed for (U-Th)/He dating. Rocks were crushed and apatite and zircon minerals were purified after magnetic and heavy-liquid separation. Apatite and zircon grains with euhedral morphology and no visible inclusions were selected under a polarized microscope. Among all these samples, only two (TGY-3 and 3.11.10.10) yielded suitable apatite grains for (U-Th)/He dating. All the samples yielded at least three single zircon grains suitable for (U-Th)/He dating. Measurements of U, Th, and He were performed at (U-Th)/He Chronology Laboratory, Institute of Geology, China Earthquake Administration. We obtained two apatite (U-Th)/He (AHe) mean ages and eighteen zircon (U-Th)/He (ZHe) mean ages (see online Supplementary Table S1 in the online Supplementary Material 1). The two AHe ages are 4.73 ± 0.73 Ma and 0.66 ± 0.07 Ma, respectively. Of the eighteen ZHe ages, fifteen are <3 Ma among which the nine ZHe ages from the vertical transect vary between a small range of ~1.0 and 1.5 Ma. The remaining three ZHe ages are >17 Ma with the oldest age of 71 Ma. Both AHe and ZHe ages are older downstream compared to ages upstream. Our ages fit into the regional pattern along this section of the Gez River where young biotite ⁴⁰Ar/³⁹Ar (ArB), muscovite ⁴⁰Ar/³⁹Ar (ArM), apatite fission-track (AFT), and ZFT ages of <4 Ma are localized in the Kongur Shan dome and older ages of >10 Ma in the lower course of the Gez River [8, 28, 31] (Figure 2).

In order to investigate the relation between the uplift of the Kongur Shan and the rang-parallel faults, we estimate an uplift model from the long river profile of the Gez River and cooling ages of multiple low-temperature thermochronometers. Our analysis is based on the detachment-limited stream power model [45, 46], which is suitable for the bedrock rivers in the Kongur Shan. In this empirically justified model, the erosion rate at any given point on the channel scales with the upstream drainage area and the local channel slope, and the elevation change at this point can be expressed as

in which $z$ is the elevation, $t$ is the time, $x$ is the distance from the outlet, $U$ is the rock uplift rate, $K$ is the erosion efficiency parameter, $A$ is the upstream drainage area as a proxy for the discharge, and $m$ and $n$ are the area and slope exponents, respectively. Given the very young cooling ages (mostly <4 Ma except near the mountain foreland) and thus rapid exhumation, the erosion rate in the Kongur Shan most likely balances against the rock uplift rate, that is, the current shape of the river’s long profile is at a steady state and experiences no changes in elevation. Thus, the stream-power equation is simplified as

Therefore, the uplift rate $U$ can be estimated for any given values of parameters $K$⁠, $m$⁠, and $n$⁠. Here, we use $m=0.45$ and $n=1$ [47-49] and the value of $K$ will be determined by inverse modeling of the thermochronological data (see below). As the uplift of Kongur Shan is likely controlled by activities on several range parallel faults (the Kongur Shan fault, Gez fault, Oytagh fault, and the MPT from west to east, Figure 2), we consider a block uplift model for which $U$ is unique within each fault-bound block but varies across the main faults, that is, the Oytagh and Gez faults.

Thermochronological ages from the Gez River are modeled with a 1D thermal model to provide constraints of the exhumation history for a potentially longer period. To estimate uplift history across the faults, our inverse modeling seeks to optimize the onset time of the current phase of uplift within the blocks between faults, that is, $t01$ for the eastern block, $t02$ for the central block, and $t03$ for the western block, as well as earlier uplift rates, $u01$⁠, $u02$⁠, and $u03$ of the eastern, central, and western blocks, respectively. We assume that the river was also in steady state prior to the current uplift phase and estimate the river elevation using

which is transformed from equation (2); $U(x)$ has values of $u01$⁠, $u02$⁠, or $u03$ according to the location along the river profile, and $z0$ is the elevation at the river outlet which is assumed to be constant at the base level.

The uplift, erosion, and topographic models described above are used to compute 1D thermal models using a finite-element method. The Gez River longitudinal profile is extracted from the SRTM DEM data with a resolution of ~30 m. The uplift is only inverted for the Gez River segment in the footwall of the Kongur Shan normal fault where river incision occurs on bedrocks. Glacier or landslide could contribute to regional denudation. However, a current V-shaped gorge implies a dominant force of fluvial erosion along this segment of Gez River. Upstream of this segment is located in the Muji-Tashkurgan Valley occupied by sediments and changes in river bed elevation should be affected by sedimentation besides river erosion. We selected the elevation of the outlet at the Tarim Basin (see Figure 2(b)) as the base level, that is, $z0$ for equation (3). The inversion process also searches for the temperature at the base of the thermal model (⁠$Tb$⁠). Other model parameters are listed in Table 1.

Then thermochronological ages are computed, using the diffusion parameters compiled in Reiners and Brandon [50] for the noble gas techniques, and the annealing models of Ketcham et al. [51] and Tagami et al. [52] for the apatite and ZFT methods, respectively. The predicted ages of all samples are compared to the observation using the log-likelihood function

in which $pi$ and $oi$ are the predicted and observed ages for sample $i$⁠, respectively, $σi$ is the $1σ$ error, and $na$ is the number of observed ages. Observed thermochronometric ages located along the river bedrock of the Gez River segment in the footwall of the Kongur Shan fault, containing our new data and published ages of different thermochronometric systems, are used for inversion (Figure 2(b); online Supplementary Tables S1 and S2). Only thermochronometric ages younger than 20 Ma are used for inversion since we have poor resolution due to the limited number of ages of >20 Ma.

In the modeling, we account for the potential spatial and temporal variances in rock uplift rates. Thus, we divide the modeling domain into three blocks separated from west to east by the Gez fault and the Oytagh fault. Rock uplift rate is allowed to vary between blocks but not within them. The final best-fit scenario, based on the constructed 1D thermal model (Table 1), is determined by minimizing the differences between the observed river profile and thermochronometric ages and the predicted ones as described above. Note that most of the thermochronometric ages are localized between the Kongur Shan fault and the Gez fault in the upstream of the modeled Gez River segment and only four ages are located beyond that region. However, since the evolution of the river profile in the upstream portion and thus the rock cooling path are closely related to elevation changes in the downstream, the final best-fit model is dependent on all the thermochronometric data along the modeled river segment irrespective of their spatial distribution.

The model optimization uses a Markov Chain Monte Carlo method [53], by which best-fit values for the unknown parameters (Table 2) are estimated. The basal temperature of the thermal model and the end time of the initial uplift are searched in the linear space, whereas the erosion efficiency $K$ and erosion rates are in the logarithm space. A total of 100,000 forward models were run, with the first 50% discarded as burn-in.

The modeling results indicate that all sampled parameter values fall within relatively stable ranges in the post-burn-in ensemble (Figure 3), suggesting successful convergence to the optimal regions of the parameter space. The best-fit model reproduces the thermochronological ages reasonably well (Figure 4).

The modeled uplift rates of the three blocks (Figure 5) suggest that no significant uplift occurred in the lower part of the Gez River prior to the Pliocene, implying that the MPT as well as the Oytagh fault was inactive until at least 5–2 Ma. Consequently, regions east of the Gez fault had much lower topography in the late Miocene compared to present-day elevations (Figure 4(a)). In contrast, the optimized models indicate that the uplift rate of the high Kongur Shan has been rapid since at least the mid-Miocene (~8 Ma), suggesting a longer history of deformation in the footwall of the Kongur Shan normal fault.

Although our inversion does not provide a tight constraint on the onset of rock uplift between the Kongur Shan normal fault and the Gez fault, it indicates that the footwall of the Kongur Shan normal fault has experienced sustained high rock uplift since ~8 Ma with a rate of ~3 mm/yr. On the contrary, regions beyond the Kongur Shan dome to the east have experienced increased rock uplift rate no earlier than 2 Ma, with a rate of ~1.7 mm/yr in the central block and ~1 mm/yr in the eastern block (Figure 5).

Our inverted rock uplift history of the footwall of the Kongur Shan normal fault implies that this fault became active at least since 8 Ma. This is slightly older than the onset time of normal faulting constrained by detrital ZFT ages of the modern river sediments collected from the Gez River catchment (~6–4 Ma) [20]. We infer that these young detrital ages <6 Ma recorded by the modern sediments from rapidly exhumed dome area mostly reflect the currently exposed bedrock ages in the drainage basin and that the preexposure history may not be completely documented by the modern sediments as they occupied only a small proportion of the age population [20].

North of the Gez River transect in the footwall of the Kongur Shan normal fault, multi-diffusion modeling of the Ar-Ar age spectral of a bedrock sample immediately below the fault surface shows rapid cooling beginning at ~8–7 Ma [9]. Our results could be compatible with this rapid cooling age that has been interpreted to date as the initiation of normal faulting [9]. And along the western boundary of the Muztagh Ata dome, ZHe and AHe ages combined with three-dimensional thermokinematic models suggested the initiation at 6.5 Ma of the southern Kongur Shan fault [33]. Given that these time estimates are not significantly different, we infer that, if not earlier, the initiation of extension in the central portion of the Kongur Shan normal fault (~8 Ma) is broadly synchronous with that in its northern and southern portions, though it is slightly earlier than that revealed by detrital ZFT and AFT results (~6–5 Ma) [8].

Our inversion provides spatially and temporally variable rock uplift rates in the Eastern Pamir where sustained high rock uplift rates localize in the Kongur Shan dome. This rock uplift pattern contributes to the decreased elevation from Kongur Tagh toward the Tarim Basin. The highest rock uplift to the west of the Gez fault within the Kongur Shan gneiss dome defines the highest topography and causes the exposure of the gneiss dome. With the sustained high rock uplift rate in the gneiss dome, a total rock uplift of ~24 km along the river channel over the last 8 Ma can be estimated. Assuming a gentle prefaulting topography and current relief of ~4–5 km between the river channel and peaks, a minimum exhumation at the mountain peaks of ~19–20 km can be estimated. This amount of rock uplift (~24 km) is sufficient to expose the gneiss from ~650°C to 700°C in the middle crust at ~9 Ma to the surface now at ~3300 m [9] if a geothermal gradient of ~30°C/km is assumed. This inferred geothermal gradient is in fact close to the value estimated by our inversion: the best-fit basal temperature of 687°C (Table 2) and model thickness of ~23 km (thickness of 20 km from base to the sea level [Table 1] plus the current elevation of ~3300 m) generate a geothermal gradient of ~30°C/km. Our thermochronometer-and-topography-derived rock uplift rate (~3 mm/yr) of the Kongur Shan dome is also in good agreement with previous exhumation rates [28, 29], revealing a rapid and protracted uplift history along with strong exhumation in the Eastern Pamir gneiss dome since ~8 Ma. The Kongur Shan normal fault is still active currently and a vertical rate of ~2.2 mm/yr [54] has been reported along the fault to the north of the Kongur Shan dome and a rate of ~1–2 mm/yr to the south of the Muztagh Ata dome [41], although both rates are lower than our inverted long-term uplift rate.

Uplift beyond the Kongur Shan dome to the east may be related to the thrusting of the Oytagh fault and the MPT. The MPT accommodated convergence between the Pamir and the Tarim during the Miocene time [43, 55], while the Oytagh fault, which merges with the MPT to the north and Kashgar-Yecheng transfer system (KYTS) to the south [8], might trigger the thrusting-related cooling in the late Miocene [20]. However, both faults do not show clear evidence of movement in the last 5–2 Ma. For example, deformation in the northeastern Pamir margin propagated northward since ~5 Ma [56], featured by the initiation of activation of the Pamir Frontal Thrust with piggy basins formation in its hanging wall at ~6–5 Ma [56-58]; the Kusilaf fault within the KYTS, which connects to the Oytagh fault, has been inactive since at least 5–3 Ma [30]. Accordingly, the MPT and the Oytagh fault might have ceased since ~5 Ma, and the current GPS velocity across the MPT shows no obvious step [18]. Thus, we infer that the increased uplift rate in the central and eastern blocks since the Pliocene is not accommodated through thrusting along single faults. It is rather associated with regional tectonics.

The late Miocene initiation of uplift and normal faulting in the Eastern Pamir postdates the initiation of dome exhumation in the Central and South Pamir [10, 11, 13]. Previous research suggested that crust thickening of the Pamir happened during ~37–20 Ma due to the India-Eurasia collision documented by prograde metamorphism in the Central and South Pamir gneiss domes as well as the Muztagh Ata dome in the Eastern Pamir [11, 59-62]. After the thickening of the Pamir crust, the transition to exhumation of these gneiss domes started at ~20 Ma [6, 11, 13, 62, 63] related to gravitational collapse of the thickened Pamir crust [11]. Shortening resumed at ≥12 Ma due to the continuous northward propagation of the Indian lithosphere [13]. Crustal thickening at this stage may have triggered peak metamorphism in the Kongur Shan dome around ~9 Ma [9]. While the axes of the gneiss domes in the Central and South Pamir are oriented perpendicular to the direction of shortening, the axis of the Kongur Shan dome is aligned with it. This geometry suggests a distinct mechanism for the formation and exhumation of the Kongur Shan dome. The localized high rock uplift observed in the Kongur Shan dome may be explained by the presence of a flat–ramp–flat thrust fault at depth in the Eastern Pamir [33], which could account for its formation, uplift, exhumation, and the subsequent eastward propagation of deformation. In this tectonic setting, the Kongur Shan dome is positioned above the ramp, which promotes doming, sustained uplift, and exhumation. This configuration also permits the eastward propagation of deformation along the flat toward the Tarim Basin. Therefore, the Kongur Shan normal fault may result from the upward extrusion of material above the ramp.

In this study, we report new thermochronometric ages across the Eastern Pamir along the Gez River. These ages, combined with existing data, were inverted together with the river’s longitudinal profile to reconstruct the uplift history. Our inversion provides spatiotemporally variable rock uplift rates across the Eastern Pamir. Sustained high rock uplift rates (~3 mm/yr) localize in the Kongur Shan dome which is due to the onset of the Kongur Shan normal fault at ~8 Ma. To the east of the dome structure, no significant uplift is observed until an increase of uplift around the Pliocene. Such an increased uplift rate is not accommodated through thrusting along single faults, given no clear evidence of motion along the MPT or the Oytagh fault. Considering the uplift history and its spatial patterns in the Eastern Pamir, we interpret that the Kongur Shan normal fault may be a result of the upward extrusion of crust materials along a flat–ramp–flat thrust fault at depth.

Apatite and zircon (U-Th)/He ages supporting this work are available in the Supplementary Material and via Zenodo: https://doi.org/10.5281/zenodo.7638540. Topographic data used in this study are available from www.hydrosheds.org. Topo-ToolBox 2 [65] was used for processing topographic data. MATLAB implementations of the affine invariant ensemble sampler (github.com/grinsted/gwmcmc) were used for MCMC sampling.

The authors declare that there is no conflict of interest regarding the publication of this article.

This project was supported by the National Natural Science Foundation of China (41961134031, 41972214).

We thank the reviewers for their constructive comments and suggestions which helped improve the article significantly.

The supplement material provides details of new data reported in this study including sample locations and single apatite and zircon (U-Th)/He ages (Table S1), and a list of published ages used for inversion in this study (Table S2).