Tectonics drives rapid exhumation of the western Himalayan syntaxis: Evidence from low-temperature thermochronometry of the Neelum valley region, Pakistan

The western syntaxis of the India-Asia collision zone is arguably the most active and spectacular example of a collisional orogeny. It features extreme rates of both mountain building and erosion. Nanga Parbat, the ninth highest mountain on Earth, is at its center, and one of the largest drainage systems in Asia, the Indus River, crosses the region (Fig. 1A). The syntaxis is surrounded by several crustal-scale thrusts (Main Mantle thrust, MMT; Main Central thrust, MCT; Main Boundary thrust, MBT) (Burg et al., 1996; Guillot et al., 2008; Wilke et al., 2012; Pêcher et al., 2002, 2008; Zeitler et al., 2001a, 2001b; Zeitler, 1985; Fig. 1B). The highest rock uplift and erosion rates have been recorded here (Zeitler, 1985; Zeitler et al., 1982; Wilke et al., 2012) with fluvial incision rates as high as 12 mm/yr (Burbank et al., 1996). The core of the syntaxis records unusually young high-temperature cooling ages and contains leucogranites that are indicative of the extreme thermal events in the region (see review of Schneider et al., 2001; Koons et al., 2013).

The recognition of these extreme rates of mountain building and destruction has sparked an active discussion regarding their causes. Two opposing models have been developed: one hypothesis considers that both rock and surface uplift are directly related to the geodetic rates of contraction (>10–15 mm/yr) (Khan et al., 2008; Seeber and Pêcher, 1998; Garcia-Castellanos and Jiménez-Munt, 2015), the high erosion rates being the consequence. The alternative posits that erosion by rivers provides extreme denudational unloading and that the high uplift rates are the consequence, implying a direct feedback between deep-seated tectonic and surface processes (Koons et al., 2013; Zeitler et al., 2001a, 2001b). This recognition has led to the formulation of the tectonic aneurysm model (e.g., Zeitler et al., 2001a, 2001b), which implies localized surface uplift and exhumation in the syntaxis region only. While this model has found wide appeal for the interpretation of the western syntaxis, studies in the eastern syntaxis question the universal application of this model (King et al., 2016; Wang et al., 2014). Recent studies have shown that in addition to the normal convergence between India and Asia the extreme topography is sustained by the horizontal mass flux into the region that occurs due to the unusual geometry of the syntaxis (Whipp et al., 2014).

Many recent studies of the western syntaxis region discuss models that explain the inferred aneurysm (e.g., Koons et al., 2013). However, the southern limits of the rapid exhumation (or the aneurysm) remain poorly constrained (Fig. 1C) because data from the critical region directly south of the Nanga Parbat syntaxis (around the Neelum valley region of Azad Jammu and Kashmir, Pakistan) are missing (Figs. 1C and 2). Data from this region are crucial if we are to discern if the region of rapid exhumation is spatially confined to Nanga Parbat region (i.e., it forms an isolated aneurysm), or if it forms a lobe around the syntaxis, thus paralleling the major structures.

In order to understand the exhumation processes at the western syntaxis and to identify the southern limits of the rapidly exhuming region, we have performed a study of the low-temperature cooling history (e.g., Reiners and Brandon, 2006) of the Neelum valley region. We report apatite U-Pb thermochronometry (recording cooling through ∼450 °C; Chamberlain and Bowring, 2000), apatite fission track (AFT) ages (recording cooling through 80–140 °C), as well as apatite and zircon (U-Th-[Sm])/He ages (recording cooling through 40–80 °C and 140–205 °C intervals, respectively), for samples from two altitudinal profiles and an along-river transect. The data are complemented by geomorphic analyses of selected fluvial channels in the Neelum River catchment in order to constrain the geomorphological and geodynamic evolution of the Nanga Parbat syntaxis.

The Himalayan mountain belt formed due to the collision of the Indian and Asian plates that initiated in the Paleocene (Ding et al., 2016). In the northwest Himalayas the main rock units are separated by regional-scale thrusts (Fig. 1B). North of Nanga Parbat, the MMT demarcates the southern boundary of Asia and the Kohistan-Ladakh arc and outlines the Nanga Parbat syntaxis. The High Himalayan Crystalline (HHC) rocks in the footwall are mainly high-grade metamorphic rocks intruded by granites (Calkins et al., 1975; Fontan et al., 2000). The HHC rocks were thrust onto the Lesser Himalayan rocks along the MCT during the early Miocene (Hodges et al., 1992). The Lesser Himalayan units are composed primarily of low-grade metamorphic and sedimentary rocks. Further south, the MBT forms the Hazara-Kashmir syntaxis (HKS), where rocks of the Lesser Himalayas are thrust onto the Neogene Sub-Himalayan molasse sediments (Kazmi and Jan, 1997; Shah, 2009).

Deformation has generally propagated from north to south. Asia and the Kohistan island arc were sutured in northern Pakistan between 102 and 85 Ma along the Main Karakoram thrust (Treloar et al., 1989). Further south, the collision of India with Asia and the Kohistan-Ladakh island arc was not synchronous: India collided in northern Pakistan with the Kohistan-Ladakh island arc along the MMT at 65–55 Ma (Kazmi and Jan, 1997), but only at 54 Ma in the Ladakh region of northwest India (Patriat and Achache, 1984). However, more recent studies have shown that initial collision occurred between India and the Kohistan-Ladakh island arc ca. 50 Ma followed by collision of this assembled landmass with Asia ca. 40 Ma (Bouilhol et al., 2013). Regardless of the chronology, this resulted in sandwiching of the Kohistan-Ladakh island arc between Asian (Karakoram) and Indian (Higher Himalayan) crust at the northwestern edge of Indian plate. Elsewhere, Indian plate is sutured directly to the Asian (Eurasian) plate along the Tsangpo suture (Bouilhol et al., 2013). HHC rocks were exhumed from mid-crustal levels along the MCT between 24 and 10 Ma (Searle et al., 2008). Deformation then shifted to lower structural levels to the south. Movement on the MBT initiated in the late Miocene (before10 Ma) in the western Himalaya (Meigs et al., 1995), followed by the development of the Main Frontal thrust (or Salt Range thrust), the youngest of all the regional faults, at 5.4–1.0 Ma (Kazmi and Jan, 1997).

Despite this well-defined sequence of major tectonic activity along the major thrusts, it is not clear if the final exhumation of rocks and the elevated topography of the present-day landscape are related to the same series of thrusting events. For example, although the main movement along the MCT is thought to have taken place during the early to middle Miocene (Butler et al., 1989; Vannay et al., 2004), AFT ages from HHC rocks in the syntaxis region are substantially younger, and an isolated area of extreme exhumation has been recognized (Zeitler et al., 1982; Zeitler, 1985). Generally AFT cooling ages are ca. 1 Ma in the core region of the Nanga Parbat syntaxis, whereas they are as old as 6 Ma toward the southwest in the Kaghan and Indus valleys (Fig. 1C). Wilke et al. (2012) reported AFT ages of ca. 10 Ma or younger from the upper Kaghan valley, but their samples are not from valley bottoms so they do not record the most recent part of the thermal history. Although the MMT may only have remained active from early Eocene until late Oligocene or early Miocene (DiPietro and Pogue, 2004), many of the low-temperature cooling ages from the Kohistan arc are significantly younger (20–10 Ma; e.g., Zeitler, 1985). Thus, much of the final stage of exhumation occurred after the majority of the tectonic history identified by geochronology, and paleomagnetic, metamorphic, and stratigraphic methods (Treloar et al., 1989; Kazmi and Jan, 1997; Patriat and Achache, 1984; Searle et al., 2008; Meigs et al., 1995).

It has been suggested that the rapid exhumation may be related to rock uplift over a crustal-scale ramp in the basal décollement, and not to individual structures (Bollinger et al., 2004, 2006), and the pronounced three-dimensional (3D) nature of the exhumation history has been documented (Whipp et al., 2007; Herman et al., 2010). While rock uplift over an unexposed basal décollement implies a more distributed young exhumation and a complicated relationship between cooling and exhumation, surface uplift of the region appears to have occurred by the Oligocene. The Deosai Plateau to the northeast of the study region is considered to be a remnant of this early surface uplift (van der Beek et al., 2009). While the surface uplift at Nanga Parbat west of the Deosai Plateau is likely to have occurred in the past few million years (Whipp et al., 2014), the Kashmir Basin south of Nanga Parbat has remained at a low elevation, accommodating terrestrial sediment for the past 4 m.y. (Burbank and Johnson, 1982).

The geology of the Neelum valley region was mapped by Malik et al. (1996) and Fontan et al. (2000) (Fig. 2). The region is in the hanging wall of the MBT and forms three distinct tectonic units separated by two thrusts, the Gumot shear zone (GSZ), a 100–200-m-wide ductile shear zone inside the HHC, and the MCT. From the structurally highest to lowest levels the tectonic units separated by these thrusts are the Shontar unit in the hanging wall of GSZ, the Kalapani unit between the GSZ and MCT, and the Salkhala unit below the MCT (inset, Fig. 2). The position of the MCT in the Neelum valley region is subject to debate. Fontan et al. (2000) proposed that the MCT crosses the Neelum River near Doarian village (Fig. 2). In contrast, Greco and Spencer (1993) mapped a mylonite zone along the MCT that extends through the Leswa section, crossing the Neelum River between Nauseri and Jura (Fig. 2). This location is widely acknowledged to represent the position of the MCT in the Neelum valley region (Wilke et al., 2012; Pêcher et al., 2008; Searle et al., 2008; Guillot et al., 2008; DiPietro and Pogue, 2004; Greco and Spencer, 1993), and is shown in Figure 1. During field work for this study we observed a mylonite zone at a similar location along the Leswa road. The timing of thrusting on these structures is not known and is correlated with movement on the MCT in other portions of the northwest Himalayas where it was active until 16 Ma (Vannay et al., 2004), and possibly reactivated later (Catlos et al., 2001; Hodges et al., 2004). The morphology of the Neelum River catchment is similar to the Indus valley to the northwest. It spans an elevation range from 700 m (at the junction with the Jhelum River) to ∼5000 m, and crosses almost the entire morphological range of the Himalayas from the Punjab plains to the high Deosai Plateau, which is thought to be a preserved remnant of the Eocene Tibetan Plateau in the east (van der Beek et al., 2009).

Rock samples were collected along high-relief sections from near the Neelum River to the adjoining peaks in the Neelum valley. The samples were collected near major regional-scale structures (e.g., the MCT and MBT) along two altitudinal profiles (Fig. 2); one near Kel in the upper parts of the Neelum valley and the second along the Leswa road in the lower part. The Leswa road leaves the main Muzaffarabad-Kel road in the Neelum valley, crosses a pass, and joins the valley again near Maira (Fig. 2). Sampling sites cover the maximum vertical range with minimum horizontal separation between samples of the same profile. Samples were also collected between the two transects along the Neelum River.

We collected 10 samples between 2000 m and 4200 m elevation along the Kel profile. The average vertical separation between samples is 200–250 m (Table 1). Seven samples were collected along the Neelum River between the Kel and Leswa sections with ∼100 m vertical separation. Thirteen samples were collected along the Leswa altitudinal profile. Samples were taken from both sides of the Leswa pass between ∼1200 m at the Neelum River and ∼2700 m at the pass (Table 1). Two samples are sandstones of the Murree Formation, one sample is Panjal andesite, and the remaining 10 samples are metamorphic rocks from the Salkhala Group (sensu Fontan et al., 2000). The quality and quantity of apatite were variable, so only a subset of samples could be used for analysis. Where possible we have combined apatite U-Pb ages, AFT ages, and (U-Th-[Sm])/He zircon and apatite ages in order to obtain the cooling history from ∼450 to ∼40 °C.

We selected 16 samples for analysis. Apatites from 12 samples were dated using the U-Pb and AFT methods (Table 2): 4 samples from the Kel profile, 3 samples from the Neelum River, and 5 from the Leswa profile. Eight samples (largely a subset of the 12 used for U-Pb and AFT dating) were dated by the apatite U-Th-[Sm]/He (AHe) method (Farley, 2002) (Table 3) and seven samples were dated by the zircon U-Th/He (ZHe) method (Reiners, 2005) (Table 4). (For details of laboratory methods, see the GSA Data Repository Item¹.)

The thermal history of Kel and Leswa profiles have been modeled using QTQt version 5.3.3 (Gallagher, 2012). Modeled thermal histories used all thermochronometry data. This incorporated the fission track annealing model of Ketcham et al. (2007) and the effect of radiation damage on He ages using the models of Flowers et al. (2009) for apatite and Guenthner et al. (2013) for zircon. The high-temperature constraint employed the weighted mean U-Pb age (28 ± 6 Ma) of all the samples that have Cenozoic (unanchored) lower intercept ages and the corresponding U-Pb apatite closure temperature of 450 ± 75 °C. The maximum temperature change (heating or cooling rate) for all model runs was limited to 1000 °C/m.y. The present-day temperature offset between the top and bottom samples of the Kel and Leswa profiles was set at 15 ± 5 °C and 10 ± 5 °C and a present-day temperature of 5 ± 5 °C and 10 ± 10 °C was used for the top sample in each profile, respectively. The paleotemperature offset is based on the total vertical separation between the top and bottom samples of a profile section. The models were run for 200,000 iterations as burn-in (which were discarded) and the subsequent 200,000 iterations were used for the analysis.

Apatite U-Pb ages for unknowns (uncorrected for common Pb) were determined by fitting isochrons on a Tera-Wasserburg (TW) concordia to determine lower intercept ²³⁸U/²⁰⁶Pb ages (Fig. 3; see Data Repository Item). The data plot into three groups: (1) apatites that are on a well-defined discordia with Cenozoic lower intercepts that range from 32.0 ± 26.0 to 17.0 ± 10.0 Ma (2σ, unanchored) or 43.0 ± 12.0 to 18.5 ± 4.6 Ma (2σ, anchored at a ²⁰⁷Pb/²⁰⁶Pb initial value of 0.837 ± 0.005 corresponding to the Stacey and Kramers, 1975, terrestrial Pb composition ca. 28 Ma) (Fig. 3A; Table 2); (2) samples that yield Cenozoic lower intercepts ranging from 42.0 ± 11.0 to 20.0 ± 19.0 Ma (2σ, unanchored) or 48.3 ± 9.1 to 28.9 ± 4.0 Ma (2σ, anchored at 0.837 ± 0.005) but also contain older grains that can be clearly identified as a separate group (Fig. 3B). It is important to note that only the younger grains were used to fit the isochrons. The young apatite age population was also identified on the basis of the trace element composition. For example, the apatites yielding the Cenozoic lower intercept in sample N5 (Fig. 3B) are depleted in La relative to the older grains that are partly affected by Himalayan exhumation. (3) The third group is apatite that yields Proterozoic lower intercepts (Fig. 3C).

For all three groups, the first U-Pb age column in Table 2 lists unanchored lower intercept TW concordia ages where the isochron regression was constrained by the analyses alone. Lower intercept TW concordia ages are listed in Table 2 as anchored at a ²⁰⁷Pb/²⁰⁶Pb initial value of 0.837 ± 0.005, which represents the common Pb composition at 28 Ma calculated by the Stacey and Kramers (1975) crustal Pb evolution model. The weighted mean age of all the unanchored Cenozoic lower intercept ages is 27.3 ± 5.6 Ma (Fig. 3D).

AFT ages range from 7.0 ± 0.4 to 2.2 ± 0.4 Ma (1σ) (Fig. 4; Table 2). Only sample K2 had sufficient confined tracks (n = 44) to allow track length measurements. It has a mean track length of 13.36 ± 1.65 µm. The kinetic parameter, D_par, shows a range from 1.11 to 1.59 µm, while the chlorine content ranges from 0.01 to 0.74 wt%. All samples passed the chi-squared (χ²) test (Galbraith, 1981, modified for laser ablation–inductively coupled plasma–mass spectrometry fission track dating by Donelick et al., 2005).

AFT ages at the base (2015 m) and top (3916 m) of the Kel profile are ca. 4 Ma and ca. 7 Ma, respectively, and show a positive correlation with elevation. AFT ages from the Leswa profile show no resolvable relationship with elevation, ranging from 3.5 to 2.2 Ma. The ages of the six samples are indistinguishable within analytical uncertainty, indicating rapid cooling through the apatite partial annealing zone, and these ages are likely linked to rapid exhumation. Three AFT ages along the Neelum River (3.9–3.7 Ma) are indistinguishable within analytical uncertainty. Sample N8, from the base of the Leswa profile, which is also the lowest sample of the Neelum River transect, yields a slightly younger cooling age of 2.9 ± 0.3 Ma (Fig. 4).

Individual uncorrected ZHe ages range from 13.3 ± 0.3 to 5.2 ± 0.1 Ma (1σ), while corrected ages range from 20.0 ± 0.4 to 6.1 ± 0.1 Ma (1σ) (Fig. 5; Table 4). Within-sample age variation is significant, typical of many ZHe studies. This may reflect U and/or Th zonation (Dobson et al., 2008) and/or radiation damage (e.g., Guenthner et al., 2013). Average ZHe ages range from 16.5 ± 2.2 to 8.8 ± 0.7 Ma and there is no significant variation with elevation.

Single grain AHe ages range from 7.6 ± 0.2 to 1.4 ± 0.1 Ma (1σ, uncorrected) and from 8.7 ± 0.2 to 2.0 ± 0.2 Ma (1σ, Ft corrected) (Fig. 5; Table 3). Average AHe ages typically overlap or are slightly younger than the corresponding AFT ages (Fig. 5). As with the AFT data, there is a weak trend of increasing AHe ages with elevation for the Kel profile, but no significant AHe age deviation with elevation for the Leswa profile is observed (Fig. 5; Table 3). Apatites from sample L6 were small and of low quality, and the oldest AHe age from this sample was not used in the calculation of the average age (Table 3).

In order to improve upon the interpretation of the ages described here, a thermal history for the region was modeled with QTQt (Gallagher, 2012), using the combined data sets from the different thermochronometers. For the Kel profile, which represents the structurally highest rocks, five samples (K1, K2, K3, K6, K10; Tables 2–4) from which we obtained AHe, AFT, and ZHe data were used for thermal history modeling. Because the vertical separation of the Kel profile samples is more than 2000 m (Table 1), a paleotemperature offset between the top and bottom samples of 60 ± 30 °C was used in the modeling. The modeled thermal history (i.e., the weighted mean model) of the Kel profile (Fig. 6) shows that the region has a complicated cooling history. The data record the cessation of the first cooling phase ca. 19 Ma, the upper sample cooling to ∼120 °C at that time. This was followed by reheating of ∼60 °C until ca. 10 Ma, followed in turn by a discrete cooling event between 10 and 8 Ma with a cooling rate of ∼100 °C/m.y. A final phase of slow cooling to surface temperatures may have involved a period of minor reheating between 8 and 6 Ma.

For the Leswa profile near the MBT, the AHe, AFT, and ZHe data from six samples (N8, L2, L6, L7, L9, L10; Tables 2–4) were used for thermal history modeling. As this sample set spans the widely acknowledged trace of the MCT (e.g., Greco and Spencer, 1993), we modeled the samples from either side of this structure separately. This yielded results that are similar to the modeled thermal history of all the samples together and thus provides evidence that all samples across the trace of the MCT have the same thermal history below the ZHe closure temperature (∼180 °C). The ∼1350 m of vertical separation between the top and bottom samples required a 40 ± 20 °C paleotemperature offset. It is interesting that the modeling results for the Leswa profile show a thermal evolution similar to that of the Kel profile, particularly since ca. 10 Ma (the timing of initiation of the MBT). The Leswa region underwent a rapid cooling event ca. 9 Ma, during which the uppermost samples cooled from ∼180 to 60 °C in <1 m.y. (Fig. 6). Samples were reheated by ∼60 °C from ca. 8 to ca. 3 Ma, followed by a second short-duration (∼1 m.y.) rapid cooling event that brought rocks from 120–100 °C to surface temperatures ca. 2 Ma.

In order to characterize the modern erosional regime we report a morphological analysis of selected fluvial channels in the region. Channel profile analysis was carried out in ArcMap 10.1 using the 1 arc second (∼30 m) resolution ASTER GDEM version 2 (Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model; https://asterweb.jpl.nasa.gov/gdem.asp) data (Fig. 7). Following common procedures, we study channel characteristics and compare them to an assumed geomorphic steady state, where channel slope (S) and drainage area (A) are related by a power law function S = k_sA^–θ, where k_s and θ are referred to as steepness and concavity indices, respectively (e.g., Wobus et al., 2006, and references therein). The values of k_s and θ can be measured directly from digital elevation model data by performing a regression analysis on a log-log plot between A and S. The normalized steepness index (k_sn) is calculated using a fixed reference concavity value (θ_ref), which is usually ∼0.45 (Wobus et al., 2006). The k_sn is proportional to erosion rate and can be used to infer zones of variable exhumation related to fault activity (Whittaker, 2012). The k_sn of all the streams was calculated in Matlab R2011b by linking it to the ArcMap 10.1 through the Profiler toolbar using a value of 0.45 for θ_ref (for details, see Whipple et al., 2007).

The k_sn values of almost all the streams are lowest near the headwaters (k_sn < 200 m^0.9) and in the lowest reaches of the Neelum River (yellow parts of streams in Fig. 7A). In contrast, the central parts of the Neelum River are characterized by high k_sn (k_sn > 500 m^0.9; red parts of streams in Fig. 7A). The transition from high to low k_sn values along the Neelum River main channel occurs at the location where it crosses the MBT and MCT (sensu Greco and Spencer, 1993). In contrast, the Kunhar and Jhelum Rivers show comparatively low k_sn values (200–500 m^0.9) in main channels except small segments of higher k_sn values that are located primarily close to the regional faults and are thus linked to changes in lithology (Figs. 7A and 1B). In summary, the variable steepness index of the Neelum River catchment has all the characteristics of a catchment that is either not in a geomorphic steady state, or a steady-state catchment with variable rock uplift rates in different subregions.

Channel profiles of the four rivers are shown in Figure 7B. These rivers provide key information for our following interpretation. The Neelum River is of principal interest to our study. The Jhelum River, which joins the Neelum from the east just below the MBT, is the principal river draining the Kashmir Basin to the south of the study region, and was chosen because the history of its confluence will be inferred from the thermochronological data. The Kunhar River is roughly parallel to the Neelum valley to the west and north. It was chosen because its confluence with the Jhelum is only 10 km downstream of the Jhelum-Neelum confluence and because it is expected to exhibit similar morphological channel characteristics. The Sind River that flows westward through the Kashmir Basin was chosen because it does not cross any of the major faults that dissect the Neelum valley region and can thus be used for comparison.

Knickpoints have been identified on the long profiles of Kunhar, Jhelum, Neelum, and Sind Rivers (Fig. 7B). The confluence of the Neelum and Kunhar Rivers with the Jhelum River occur within ∼10 km of each other on the western limb of a major antiformal structure of northwestern Himalayas, the HKS. These three rivers show knickpoints along their long profiles at a distance of ∼90 km to ∼110 km upstream from the confluence. It is interesting that the Sind River, which does not cross any of the major Himalayan structures and joins the Jhelum River in the Kashmir Basin, also has a knickpoint ∼80 km upstream, where the river long profile shows a change from apparently equilibrated channel (downstream) to a stepped morphology (upstream). Knickpoints identified here show a good match with the k_sn values of the corresponding rivers (Fig. 7).

The thermochronometric analyses provide series of ages that document different parts of the cooling history. Starting with the highest closure temperature system (U-Pb apatite) it was shown that the samples plot into three distinct age groups, and we interpret these three groups in terms of the effects of cooling due to exhumation related to the Himalayan orogeny as (1) completely affected, (2) partially affected, and (3) unaffected ages. One example of each is shown in Figure 3 (all other data are listed in the Data Repository Item).

All samples from the Kel profile in the hanging wall of the GSZ (according to the map of Fontan et al., 2000) are within the HHC and show completely affected apatite U-Pb ages that range between 43.0 and 18.5 Ma. Assuming a surface temperature of 20 °C, a closure temperature of 450 °C for U-Pb system in apatite, and following Wilke et al. (2012) by assuming 30 °C/km as the geothermal gradient, a total exhumation of more than 14 km is inferred for the structurally highest sample and more than 16 km for the structurally lowest sample of the Kel profile. The partially affected samples of the second group all come from the lower part of the Neelum valley and consist of granites to granitic gneisses (Table 1; Figs. 2 and 3E). These analyses may represent old (Proterozoic?) apatites that have been partly to completely affected during Himalayan orogenesis, with possible additional growth of new metamorphic apatite. The third (oldest) group is from the Leswa profile, where four gneisses yield Proterozoic apatite U-Pb ages implying that they have not been affected by Himalayan tectonothermal events (Fig. 3E; Fig. S1 in the Data Repository Item). These old U-Pb apatite ages suggest that the rocks of the Lesser Himalayas (between the MBT and MCT) have undergone total Himalayan exhumation of <∼14 km but certainly >5 km (based on the late Cenozoic AHe, AFT, and ZHe age data).

The low-temperature thermochronometry data indicate that all samples record Himalayan tectonics. Both profiles record cooling through the ZHe partial retention zone (180–140 °C) between 17 and 9 Ma, with no significant difference between the two transects. Assuming generally accepted values for the closure temperatures of the different systems (i.e., ZHe = cooling through 180 °C and AHe = cooling through 50 °C), a mean cooling rate can be estimated. The average ZHe age is ca. 11 Ma and the average AHe age is ca. 4 Ma, reflecting a cooling rate for the region of ∼19 °C/m.y., broadly consistent with the thermal history model derived from QTQt. This is consistent with an exhumation rate of >0.6 mm/yr (for an assumed regional geothermal gradient of 30 °C/km; Wilke et al., 2012). However, due to the effect of heat advection that accompanies rapid exhumation, these calculated exhumation rates should be considered as an upper limit.

This new data set allows us to propose new constraints on the history of rock uplift and exhumation in the Neelum valley region and consequently resolve the evolution of the western Himalayan syntaxis region. While previous studies have developed 3D models, we show that much of the thermal history derived here can be understood using a simple 2D approach, based on a southwest-northeast profile from the Jhelum River confluence toward the Deosai Plateau (Fig. 8).

The first constraint on the tectonothermal evolution is given by apatite U-Pb ages that constrain exhumation through ∼450 °C (Chamberlain and Bowring, 2000). These data show that there are regions that have not undergone substantial Himalayan tectonothermal events, while other samples have been partly or fully affected by Himalayan tectonics in the Cenozoic (Fig. 3E; Table 2). The boundary between these two domains is between samples L2 and L6 along the Leswa road. This likely corresponds to the MCT, and suggests that movement on this structure in the Oligocene (ca. 27 Ma, as given by the weighted mean U-Pb cooling age) caused rock and surface uplift, steepening of the topography, and ultimately erosion of the hanging wall. As the samples with young apatite U-Pb ages from the Kel profile are also in the hanging wall of the GSZ, this Oligocene rock uplift may have been generated in part by motion along the GSZ.

The differing U-Pb ages along the Leswa profile are consistent with the position of the MCT as suggested by Greco and Spencer (1993; Figs. 1 and 2). We suggest that the MCT remained inactive since at least the early-middle Miocene, as thermal history modeling of samples from either side of the MCT suggests a consistent shared cooling evolution since that time. The mixed U-Pb age samples record the combined effects of thrusting, initially along the MCT and later along the MBT (Fig. 8). The young U-Pb apatite age samples from the Kel profile are the result of greater amounts of exhumation and the combined effect of thrusting along the GSZ, MCT, and MBT (Fig. 8). Detrital zircon U-Pb ages also confirmed enhanced Himalayan exhumation at 35–23 Ma (Ding et al., 2016).

Following this early stage of exhumation, the Kel profile rocks were at 180–100 °C ca. 19 Ma. Cooling to this point was likely a result of continuous erosion following movement on the MCT and GSZ so that the relief created by the earlier event was removed (Fig. 8). It is possible that MCT was active during this entire period. Reheating (∼60 °C) between 19 and 10 Ma has no obvious associated regional tectonic event; it may simply relate to the rapidly exhuming hot rocks from great depths from the south focused into the Nanga Parbat region (Koons et al., 2002) and the resultant active hydrothermal system (Chamberlain et al., 2002; Craw et al., 1994). Rapid cooling at 10–8 Ma coincided with the first movement of the MBT (Meigs et al., 1995), and likely records the onset of rock uplift and erosion of the hanging wall. This is evident from the thermal history of the Leswa profile (Fig. 6). However, rapid cooling of the Kel profile at the same time (i.e., 10–8 Ma) may represent activity along the GSZ. The north to south propagation of deformation in the Himalayas at that time (Vannay et al., 2004) is consistent with the MBT being the main active tectonic structure, which in turn implies that the MCT had ceased to be active by that time (Vannay et al., 2004). Both sides of the MCT have a similar cooling history from between 27 Ma (based on U-Pb apatite data) and 11 Ma (based on ZHe data) onward. Treloar et al. (1992) suggested that the southwest-verging Kashmir thrust system (which includes both the MCT and MBT) at the eastern limb of the HKS (southwestern part of our study region) was actively propagating to the southwest with its frontal thrusts cutting up through the south-southeast–verging Pakistan thrust system on the western limb of the HKS. This also supports the hypothesis that the MBT was the locus of the main tectonic activity ca. 10 Ma. Formation of the MBT was followed by nearly isothermal to slow reheating conditions (Fig. 6). Final cooling recorded by the Kel profile is likely due to incision of the Neelum River catchment into an existing topographic plateau, whereas the Leswa profile shows a second rapid cooling event ca. 3 Ma (Fig. 6), possibly due to thrusting along the MBT.

The interpretation of the final stages of exhumation is supported by aspects of the geomorphic analysis. For example, the high k_sn values in the hanging wall of the MBT in the Neelum valley may indicate recent motion along this structure, but the sharp lithological contrast across its trace may also explain the change in the steepness index (Figs. 2 and 7). We interpret the low k_sn parts in the headwaters of almost all the small streams as areas where the incision has not yet propagated into an uplifted region. This is consistent with the Deosai Plateau being the remnant of an Eocene rock and surface uplift event. The high k_sn in the main channel suggests that the entire Neelum River catchment may be the product of a young river piracy event by the Jhelum River, possibly due to the formation of the Kashmir Basin ca. 5 Ma (Alam et al., 2015) and the resultant drainage convergence into the Jhelum River. Rapid cooling ca. 3 Ma as documented predominantly for the Leswa profile may be related to the incision of ∼1000 m (Fig. 7B) due to the passing of a knickpoint along the Neelum River. Our interpretation is also consistent with the slight increase in AFT age toward the headwaters of the catchment.

The AFT data can be combined with previously published ages from the Nanga Parbat region to make a substantial improvement on the interpretation of the nature of exhumation around the western syntaxis. According to the interpretation of Zeitler (1985), contours of AFT age are continuous across the syntaxis and are perpendicular to the major structures in the region, with ages at base level in the Neelum valley region predicted to be older than 5 Ma (Fig. 1C). For tens of kilometers around Nanga Parbat, the topology of these contours remains consistent with the AFT data from Treloar et al. (2000), van der Beek et al. (2009), Wilke et al. (2012), and other cooling age data (e.g., Schneider et al., 2001). Zeitler (1985) inferred a region of extreme exhumation only in the region at, and north of, Nanga Parbat, which has remained the basis for many crustal exhumation models since (e.g., Zeitler et al., 2001a, 2001b; Koons et al., 2013).

In order to broaden the interpretation of the published AFT cooling ages to the larger syntaxis region yet remain compatible with previous studies, we follow the method of Zeitler (1985) and show only the low-altitude samples K10, N5, N7, and N8 (yellow dots in Fig. 9; Table 2). AFT ages of these samples range from 3.9 ± 0.4 to 2.9 ± 0.3 Ma, but given that samples from higher elevations are not substantially older (Fig. 4), we consider the new contours, and thus the interpretation, to be robust. Combining our new data with published AFT cooling ages shows that the age contours south and east of Nanga Parbat region are broadly parallel to the major thrust structures around the syntaxis. It appears that the contours west of the syntaxis likely make a sharp bend southward at the MBT. Generally, the exhumation through the apatite partial annealing zone temperatures appear to be governed by movement on the major thrust zones (MBT, MCT; Fig. 9) with the youngest ages being located near the MCT and MBT. The extent to which the tectonic aneurysm model of Zeitler et al. (2001a, 2001b) can be applied appears to be spatially limited. The distribution of higher temperature cooling ages and the occurrence of low-pressure granulites clearly document the unusual exhumation at Nanga Parbat (Koons et al., 2013), but this region may be small and cannot be extended to the regions of the Neelum valley and further south.

The fact that the AFT cooling age contours for all ages older than ca. 2 Ma are parallel to the major thrusts in the western syntaxis supports an exhumation model that advocates tectonics as the principal driving mechanism for at least the last stages of exhumation (Fig. 9). This is also supported by the slight trend in AFT ages from our low-altitude samples that increase up river, along with the ranges of AFT ages from the two profiles, Kel and Leswa. The Leswa profile is close to the main faults in the downstream parts of Neelum River and has the youngest AFT ages, whereas the Kel profile, which is further from the regional faults, shows somewhat older AFT ages (Table 2; Fig. 3). However, the final stage of exhumation is only related to the incision of the Neelum River catchment into a plateau formed ca. 10 Ma due to activity along the MBT (Fig. 8). This interpretation is supported by stream power analysis (Fig. 6), although the signal from the AHe ages is not clear. Nevertheless, it is consistent with the interpretation of van der Beek et al. (2009), who argued that the Deosai Plateau in the northeastern part of the region is a remnant of the Tibetan Plateau. As such, our interpretation may be of tectonic importance to a much larger region than investigated here.

The results from a multiple thermochronometric dating and thermal modeling approach combined with geomorphic stream power analysis allow us to conclude the following.

Apatite U-Pb data constrain the position of the MCT in the Neelum valley region, which has been the subject of controversy. Samples with Proterozoic U-Pb ages all are part of the Leswa profile. The boundary between the samples unaffected by Himalayan tectonics and all other samples is interpreted to represent the trace of MCT along the Leswa profile.

Average ZHe, AFT, and average AHe ages are generally 16–9 Ma, 7–3 Ma, and 6–3 Ma, respectively (Fig. 4). The ages from the two apatite low-temperature thermochronology systems are largely indistinguishable for both profiles, indicating exhumation rates of 1 mm/yr to 0.3 mm/yr, with only a short residence time in the apatite partial annealing zone and AHe partial retention zone.

Thermal history modeling of our data is consistent with thrusting on the MBT ca. 10 Ma in both the Kel and Leswa profiles (Fig. 7). The MCT is also suggested to have been inactive since the early-middle Miocene, as samples across the MCT trace show the same thermal history for the last 11–10 m.y., even when modeled as separate groups of hanging-wall and footwall samples.

The Neelum River is characterized by a steepness index of k_sn > 500 m^0.9, but both the lowest reaches (where it joins the Jhelum River) as well as the headwaters of the catchment have significantly lower steepness indices, k_sn = ∼20–500 m^0.9. We interpret this as evidence for strong geomorphic disequilibrium caused by recent motion on the MBT and/or MCT and continuous catchment capture in the headwaters. Knickpoints have been identified along the channel profiles of the Kunhar and Neelum Rivers at a distance of ∼90 and ∼110 km, respectively, upstream of their junction with Jhelum River (Fig. 6). This indicates a wave of erosion propagating upstream that may have started as a result of base-level fall at 5–3 Ma due to the formation of Kashmir Basin and resultant drainage convergence into Jhelum River.

In combination with AFT ages from the literature, our data define AFT cooling age contours around the entire western syntaxis region. The contours do not close in the south of Nanga Parbat, but extend toward the southeast, trending parallel to regional-scale thrust faults (Fig. 9). The distribution of exhumation around the syntaxis thus follows the well-established pattern known from much of the Himalayan front, rather than representing an isolated region of extreme exhumation.

This work was supported by the Higher Education Commission (HEC) of Pakistan and Austrian Agency for International Cooperation in Education and Research (OeAD-GmbH) at different levels. Partial support for field work was provided by the Heinrich-Jörg Stiftung funding scheme of the University of Graz, Austria. We thank Luigia Di Nicola, Mahmoud Hassan, Gustav Hanke, Tamer Abu-Alam, Jamil Hassan, Daniel Döpke, and Katarzyna Luszczak for providing help at various stages. David Wallis, Peter van der Beek, and four anonymous reviewers helped to improve earlier versions of this manuscript. We thank Damian Nance for his editorial handling of the manuscript.