Abstract

Competing hypotheses suggest that Himalayan topography is sustained and the plate convergence is accommodated either solely along the basal décollement, the Main Himalayan thrust (MHT), or more broadly, across multiple thrust faults. In the past, structural, geomorphic, and geodetic data of the Nepalese Himalaya have been used to constrain the geometry of the MHT and its shallow frontal thrust fault, known as Main Frontal thrust (MFT). The MHT flattens at depth and connects to a hinterland mid-crustal, steeper thrust ramp, located ∼100 km north of the deformation front. There, the present-day convergence across the Himalaya is mostly accommodated by slip along the MFT. Despite a general agreement that in Nepal most of the shortening is accommodated along the MHT, some researchers have suggested the occurrence of persistent out-of-sequence shortening on interior faults near the Main Central thrust (MCT).

Along the northwest Himalaya, in contrast, some of these characteristics of central Nepal are missing, suggesting along-strike variation of wedge deformation and MHT fault geometry. Here we present new field observations and seven zircon (U-Th)/He (ZHe) cooling ages combined with existing low-temperature data sets. In agreement with our previous findings, we suggest that the transect of cooling age patterns across the frontal Dhauladhar Range reveals that the Main Boundary thrust (MBT) is a primary fault, which has uplifted and sustained this spectacular mountain front since at least the late Miocene. Our results suggest that the MBT forms an ∼40-km-long fault ramp before it soles into the MHT, and motion along it has exhumed rocks from depth of ∼8–10 km. New three-dimensional thermokinematic modeling (using Pecube finite-element code) reveals that the observed ZHe and apatite fission track cooling ages can only be explained by sustained mean MBT slip rates between ∼2.6 and 3.5 mm a–1 since at least 8 Ma, which corresponds to a horizontal shortening rate of ∼1.7–2.4 mm a–1. We propose that the MBT is active today, despite a lack of definitive field or seismogenic evidence, and continues to accommodate crustal shorting by out-of-sequence faulting. Assuming that present-day geodetic shorting rates (∼14 ± 2 mm a–1) across the northwest Himalaya have been sustained over geologic time scales, this implies that the MBT accommodated ∼15% of the total Himalayan convergence since its onset. Furthermore, our modeling results imply that the MHT is missing a hinterland mid-crustal ramp further north.

INTRODUCTION

Various orogenic growth models have been proposed to explain the evolution of the central Himalayan tectonic system, spanning more than 1000 km from the Sutlej-Beas region of northwest India into central Nepal, since the Miocene. In general this model of central Himalaya has been used to explain how the Himalaya deform; it is characterized by the exposure of the Himalayan high-grade metamorphic core, which has been extruded along subparallel shear zones, the Main Central thrust (MCT) and South Tibetan detachment, active between the early and middle Miocene. This unit terminates in the northwest Himalaya (e.g., DiPietro and Pogue, 2004; Vannay et al., 2004; Yin, 2006; Webb et al., 2011, and references therein). The central Himalaya (Fig. 1A, inset) is topographically characterized by a sharp physiographic transition (PT2) (e.g., Hodges et al., 2001; Wobus et al., 2003) between the Lesser Himalayan foothills to the south, having moderate relief and topographic peaks mostly ≤3000 m, and the steeply rising High Himalaya to the north, having high relief and peaks exceeding 6000 m (Lave and Avouac, 2001). Structural, petrologic, and thermochronologic studies indicate that above the PT2, exposed high-grade metamorphic rocks have been rapidly exhumed from mid-crustal depth since the middle-late Miocene. The PT2 correlates with major knick zones in orogen-crossing longitudinal river profiles (e.g., Seeber and Gornitz, 1983; Duncan et al., 2003), high fluvial incision rates (Lave and Avouac, 2001), and a long-standing (106 yr) erosional gradient (Finlayson and Montgomery, 2003; Olen et al., 2015). This is supported by the fact that the PT2 corresponds to the transition from relatively old mineral cooling ages (zircon (U-Th)/He [ZHe] older than 10 Ma, apatite fission track [AFT] older than 5 Ma) in the foothills to young (ZHe younger than 10 Ma, AFT younger than 5 Ma) ages in the High Himalaya (e.g., Blythe et al., 2007; Robert et al., 2009; Herman et al., 2010; Thiede and Ehlers, 2013). Based on these characteristics of the central Himalaya, a debate continues regarding how the wedge is deforming. One of the key questions is why the steep and rapidly eroding front of the High Himalaya is located ∼100 km north from the southern tip of the wedge, defined by the Main Frontal thrust (MFT). This ∼30–50-km-wide zone of deep exhumation along High Himalayan front correlates with a belt of clustered microseismicity at depths of 15–25 km (Pandey et al., 1995; Mahesh et al., 2015; Elliott et al., 2016). Some have attributed rapid exhumation to a mid-crustal ramp along the main plate boundary fault, the Main Himalayan thrust (MHT) (Pandey et al., 1995; Lavé and Avouac, 2001), or underplating and duplex formation in the Lesser Himalaya (Bollinger et al., 2004; Herman et al., 2010). Conversely, others have argued for active out-of-sequence thrusting along steep splay faults maintaining the PT2 front of the High Himalaya (Wobus et al., 2003, 2005; Bai et al., 2016; Whipple et al., 2016). In the central Himalaya, the foothills south of PT2 compose a ∼50-km-wide belt of openly folded Lesser Himalayan rocks in the hanging wall of the Main Boundary thrust (MBT), where tectonic klippes of High Himalayan crystalline complex (HHC) or Greater Himalayan sequence are preserved locally (e.g., Shimla or Kathmandu klippe). Between the MBT and active fault systems bounding the Himalayan wedge to the south, the MFT, deformed and uplifted foreland sediments of the Sub-Himalaya are limited to a narrow (∼20 km) band. Competing hypotheses suggest that plate convergence is accommodated either predominantly on the MHT-MFT or more broadly across multiple out-of-sequence faults within the Sub-Himalaya, or MBT thrust fault arrays within the wedge in this portion of the orogen.

Domains that structurally deviate from this central Himalayan geometry call into question the application of these models to the entire mountain belt. The northwest Himalaya west from 77°E differ structurally and physiographically from the central Himalaya (DiPietro and Pogue, 2004; Yin, 2006; Deeken et al., 2011; Webb et al., 2011; Webb, 2013). Here the wedge is characterized by a strongly undulating (in map view) topographic front forming large reentrant structures bounded by the MBT, such as the Kangra and Dehradun reentrants (see Fig. 1). In the reentrants, the Sub-Himalayan domain is broadly exposed and as much as 100 km wide, whereas the exposure of Lesser Himalayan rock is limited to a width of less than 10 km. In contrast to the MBT, the MFT trace is more regular and remains arc parallel. Recent studies indicate that (1) significant crustal shortening is accommodated by out-of-sequence deformation within the Sub-Himalaya (Powers et al., 1998; Thakur et al., 2014; Dey et al., 2016), and (2) ongoing exhumation occurs in the hanging wall of the MBT (Deeken et al., 2011). Improved knowledge of the main fault geometries and activity of the wedge is essential to understand the lateral variation in growth of the Himalayan orogenic belt and where crustal shortening is accommodated within the wedge other than along the MFT.

In this study we focus on a transect across the northwestern Himalaya, 100 km northwest of the Sutlej River, the western termination of central Himalayan tectonic regime (Figs. 1B and 2; DiPietro and Pogue, 2004; Vannay et al., 2004). Specifically, we focus on the Kangra reentrant, the Dhauladhar Range, located in the MBT hanging wall, and the Chamba region in the western Himachal Pradesh, India (Fig. 1). We report seven new ZHe ages that complement previously published ZHe and AFT ages (Deeken et al., 2011) to better constrain the temporal and spatial evolution of the 100-km-long segment within the wedge. By combining physiographic analyses, new and available thermochronologic data, and three-dimensional (3D) thermokinematic modeling, we test the applicability of varying fault geometries to explain the obtained cooling age pattern. More specifically, we constrain the geometry and depth and length of the MBT fault ramp and the MHT. Using Pecube (finite-element code) thermokinematic modeling (Braun, 2002b; Robert et al., 2011; Braun et al., 2012), we obtain better control on estimates of average MBT fault displacement rates, indicating Quaternary faulting. We constrain the spatial distribution of crustal shorting across the northwest Himalaya since the late Miocene, which differs significantly from shortening patterns in the central Himalaya.

GEOLOGIC SETTING OF THE NORTHWEST AND CENTRAL HIMALAYA

The Himalaya range is an orogenic wedge formed by stacking thrust sheets scraped off the Indian crust as it was driven beneath the Tibetan Plateau to the north. All thrust faults within the wedge are assumed to sole into the plate bounding main basal décollement, the MHT (e.g., Hauck et al., 1998; Stevens and Avouac, 2015; Gao et al., 2016). Approximately 20 mm a–1 of Indian-Eurasian convergence is accommodated across the central Himalayan wedge and ∼14 mm a–1 across the northwestern Himalaya (Kundu et al., 2014; Stevens and Avouac, 2015). The central Himalayan wedge is divided into fault-bound tectonostratigraphic units exposed from south to north. The MFT delineates the southernmost front of the wedge (Fig. 1C), forming the present-day surface-rupturing ramp of the MHT, and has accommodated the majority of crustal shortening over the past ∼2 m.y. (e.g., Powers et al., 1998; Lavé and Avouac, 2000; Jouanne et al., 2004). In its hanging wall, the Sub-Himalayan fold-thrust belt is bounded to the north by the MBT. The Lesser Himalayan rocks in the MBT hanging wall are unmetamorphosed to moderately metamorphosed Proterozoic to Eocene sedimentary cover of the Indian basement (e.g., Gansser, 1964; Hodges, 2000; Webb et al., 2011). Further to the north, the MCT separates the high-grade metamorphic HHC complex. In the central Himalaya, a 20–30-km-wide transition zone of microseismic activity has been recorded ∼100 km north of the surface trace of the MFT and beneath the location where the MCT is exposed at the surface (Pandey et al., 1995; Mahesh et al., 2015). In the past this has been interpreted to be located between the seismically active segment to the south and the aseismic segment to the north, and has been determined to be a ∼20° dipping hinterland mid-crustal ramp in the MHT (Cattin and Avouac, 2000; Grandin et al., 2015; Elliott et al., 2016). North of the HHC, the Sangla and Zanskar shear zones represent regional segments of the South Tibetan detachment, juxtaposing the low-grade metamorphic basal section of the Tethyan Himalayan strata, the Haimantas, over the high-grade HHC (e.g., Dezes et al., 1999).

LOCAL GEOLOGY AND FIELD OBSERVATIONS

Himalayan seismotectonic studies (Seeber et al., 1981; Rajendra Prasad et al., 2011) and balanced cross sections across the northwestern Sub-Himalaya (Yeats et al., 1992) account for what is known about the MHT fault geometry in the northwest Himalaya. The MHT fault includes several flat and ramp segments. The southernmost segment of the MHT, referred to as MFT, dips ∼15° northeast until a depth of ∼4 km, where it bends into a 3°–5° north-dipping basal décollement (Powers et al., 1998). The surface trace of the MBT at the base of the Dhauladhar Range intersects the MHT at a depth of ∼8 km below the surface. Further north, the MHT is assumed to form a flat décollement that continues 50–100 km north beneath the High Himalaya and Tethyan Himalaya (Rajendra Prasad et al., 2011). A notable scatter in microseismic activity spreads over an 80-km-wide region north of the MBT trace and may be related to the north-dipping MBT fault ramp or the MHT flat (Rajendra Prasad et al., 2011; Gahalaut and Kundu, 2012). In the Kangra sector, the Himalayan mountains rise abruptly from the Sub-Himalaya with elevations <1000 m above sea level (asl) to a series of peaks exceeding 5000 m asl forming the crest lines of Dhauladhar and Pir Panjal Ranges (Figs. 1A, 1B). Here, whereas the Sub-Himalaya is a 100-km-wide belt in the Kangra reentrant, the Lesser Himalaya units form a narrow (4–5 km) belt between the MBT and the MCT along the topographic front of the mountain range. Weakly to moderately metamorphosed sedimentary rocks of the Haimantas Group crop out in the hanging wall of the MCT in the Chamba syncline. The Haimantas Group consists of an ∼8-km-thick, continuous turbidite sequence of late Precambrian to early Paleozoic age and is correlated with the Phe Formation at the base of the Tethyan Himalaya in Zanskar (Frank et al., 1995; Fuchs and Linner, 1995, and references therein). The Haimantas Group is folded into large, open anticlines and synclines with amplitudes of several kilometers and forms the High Himalayan nappe in the MCT hanging wall (Frank et al., 1995). Exposure of high-grade metamorphic crystalline units is limited to the Kishtwar and Kullu-Rampur tectonic windows, located toward the northwest and east of our study area, respectively (Figs. 1B and 2). In contrast to the high-grade metamorphic HHC in the central Himalaya, at the Kangra reentrant the MCT juxtaposes intermediate-grade metamorphic Haimantas Group onto low-grade Lesser Himalayan sediments (Frank et al., 1995; Yin, 2006). This could be related to a lateral ramp within the MCT, such that the MCT cuts upsection from the central Himalaya to Chamba (Steck, 2003). Alternatively, the HHC in Chamba could form a southward-thinning wedge between the Tethyan Himalaya at the top and the Lesser Himalaya at the base, in which case the MCT would simply represent the original flat-over-flat relationship between the Tethyan and Lesser Himalaya when the HHC was thrust over the Lesser Himalaya (Yin, 2006).

In the Lesser Himalayan rocks exposed above the MBT (southern flanks of the Dhauladhar Range both located north of Dharamsala and Palampur; Fig. 1), we observe strongly sheared quartz veins in phyllites and prominently developed layer-parallel schistosity suggesting strongly sheared rocks at lower and moderate greenschist facies conditions (see Figs. 2B, 2C). Brittle fault planes and regional schistosity dip ∼33°–38°NE. In addition, we observed millimeter-sized garnet porphyroblasts in mica schists further up the sections and at the base of the MCT hanging wall (just a few kilometers stratigraphically upsection from the MBT; Fig. 2D), which we relate to Barrovian metamorphism prior to early Miocene MCT activity. These observations indicate an inverse metamorphic gradient across the MBT that suggests that the MBT has accommodated significant displacement.

Previous Work on Cooling History of the Dhauladhar Range and Motivation

Previous publications (Deeken et al., 2011; Adlakha et al., 2013) provide key information about the cooling history of the Dhauladhar Range; however, the sparse sampling along the base of the Dhauladhar Range front did not allow them to constrain exhumation depth. It could not be ascertained if the Pliocene ZHe ages were fully reset and if they were exhumed from below the ZHe partial retention zone, or whether they were within the ZHe partial retention zone for a prolonged period of time; the answer to this provides key information on the MBT fault geometry (length and depth) and constrains the amount of exhumation accommodated by the MBT.

METHODS

Thermochronology

We analyzed seven samples from Lesser Himalayan units, the Haimantas metasediments, and Paleozoic intrusive rocks from the MCT hanging wall along the Dhauladhar Range front (Figs. 1 and 3). Sample preparation and analytic details are provided in the Data Repository1.

3D Thermokinematic Modeling: Pecube

We have used Pecube thermokinematic modeling (Braun, 2003; Braun et al., 2012) to investigate probable time-temperature paths of the MBT hanging-wall rocks, constrain the onset of faulting, estimate mean MBT fault displacement rates, and test different geometries of the frontal segments of the MBT.

Time-temperature paths reflect the rock exhumation histories, which are controlled by the particle pathways within the orogenic wedge, i.e., they record the interaction between tectonic transport and erosion (Stuwe et al., 1994; Mancktelow and Grasemann, 1997; Braun, 2002b; Ehlers, 2005; Braun et al., 2012). For example, when rocks move over a ramp within a thrust fault, the transport vectors attain a significant vertical component, which results in rock uplift, builds topography, and drives erosion (Stüwe et al., 1994; Mancktelow and Grasemann, 1997; Braun, 2002b). The thermal field of the upper crust is therefore strongly disturbed, resulting in both vertical and lateral heat transport (Huntington et al., 2007) and can be quantified with a 3D thermokinematic model such as Pecube (Braun, 2003; Braun et al., 2012).

For Pecube modeling, we developed a series of different possible MHT-MBT-MFT geometries that satisfy the constraints of studies from Powers et al. (1998), Frank et al. (1995), Steck et al. (1993), Steck (2003), and Rajendra Prasad et al. (2011). However, especially north of the MBT trace, the MHT geometry is poorly constrained due to lack of available seismic data.

Model Setup

We used the 2012 modified version of Pecube (Braun, 2003; Braun et al., 2012), a finite element code that solves the three-dimensional heat transport equation in a crustal block affected by vertical and/or horizontal material transport. In its forward mode, Pecube predicts thermochronologic age distributions (e.g., apatite (U-Th-Sm)/He, ZHe, AFT) that result from a prescribed velocity field defined by fault geometries and slip velocities and a stable or evolving topography. Inverse modeling explores a parameter such as a range of slip velocities or fault dip angles. For each combination of parameters, the misfit between the model predictions and a cooling ages data set is computed using a chi-squared approach, which indicates the best-fitting parameter set.

The dimensions of our Pecube model are 320 × 50 × 25 km (length, width, depth); the model is oriented 25° north-northeast and includes the MBT, the MFT, and their common north-northeast–dipping low-angle sole thrust, the MHT (Fig. S1). The faults trend perpendicular to the long axis of the model, and all horizontal particle motion is parallel to the model swath.

Fault-slip rates were derived from present-day global positioning system–based convergence rates (Schiffman et al., 2013; Kundu et al., 2014; Stevens and Avouac, 2015) corresponding to 14 mm a–1 in the direction of model transport. We assume that present day convergence rates have not changed significantly since ca. 10 Ma (Molnar and Stock, 2009; Copley et al., 2010). Following previous thermokinematic modeling approaches (Bollinger et al., 2004; Robert et al., 2009; Herman et al., 2010), we split the total convergence rate of 14 mm a–1 into underthrusting of the Indian plate below the MHT and overthrusting of MBT and the MFT. In particular, we set the underthrusting rate to a constant value of 5 mm a–1 for the duration of the model runs, while the remaining 9 mm a–1 are along the MFT minus a fraction accommodated by the MBT.

We assume that the rock uplift resulting from overthrusting is compensated by denudation such that the topography remains constant throughout the model runs (Robert et al., 2009; Herman et al., 2010) and evaluate this in the model sensitivity section below. We imposed the following constraints on fault activity: (1) ZHe ages obtained from high elevations (>3000 m) across the studied range yield old ages (18–12 Ma; Figs. 1 and 3); (2) activity along the MBT was established between ca. 8 and 10 Ma (Meigs et al., 1995; Brozovic and Burbank, 2000; DeCelles et al., 2001). We assumed 10 Ma as a minimum onset time for the MBT. Model runs start at 18 Ma and are subdivided into three periods. During period 1 (18–12 Ma), the MBT is not active; all shortening is accommodated along the MFT (velocity, VMFT = 9 mm a–1). Particles that are not thermally reset during subsequent steps can yield 18 Ma cooling ages. During period 2 (12–8 Ma), the MBT starts at a velocity of VMBT = 1 mm a–1 and most overthrusting is accommodated along the MFT (VMFT = 8 mm a–1). During period 3 (8–0 Ma), the MBT full fault activity occurs at variable rates (VMBT = 1–5 mm a–1). Present-day topography was determined from the 90 m Shuttle Radar Topography Mission digital elevation model version 4.1 (Jarvis et al., 2008). We used rock thermal diffusivity measurements of Thiede et al. (2009) and inferred values for all other physical parameters of the rocks within geologically meaningful ranges (see Data Repository Table 1). The thermal diffusivity ranges between 28 and 60 km2/m.y. The mean surface temperature at elevation, Z = 0 m, was set to 20 °C, and the atmospheric lapse rate was set to 6 °C/km. (For further details, see Data Repository Figs. S2–S6).

RESULTS

Thermochronology

The new ZHe ages are obtained from elevations of 1.9–3.1 km asl and range between 3.0 ± 0.2 and 6.7 ± 1.2 Ma (Table 1). In combination with previously published ZHe and AFT cooling ages (Deeken et al., 2011), these data compose a densely sampled 80-km-long range-perpendicular transect across the hanging wall of the MBT (Fig. 1). Deeken et al. (2011) found out that AFT and ZHe cooling ages reveal marked variations in long-term exhumation rates between the humid frontal range, with very young (younger than 4 Ma), elevation-invariant AFT ages obtained in a ∼40-km-wide zone, and the semiarid orogen interior, where AFT and ZHe ages are older (11.8 ± 1.2 and 18.1 ± 1.8 Ma) and increase with sample elevation (Fig. 3). The previously published ZHe ages confirm southward-younging cooling age trends (15 ± 3 Ma in the interior to 8 ± 6 Ma in the frontal ranges). These differences have been attributed to lateral heat advection resulting from southward thrusting on the MHT and MBT faults (Deeken et al., 2011). We interpret the elevated ZHe ages (older than 10 Ma, >3.7 km asl) in the Dhauladhar and Pir Panjal Ranges (Fig. 3A) as related to exhumation and cooling due to the terminal MCT fault activity in the middle Miocene and to having stayed in the partial retention zone and above closure temperature since then. The ZHe ages, obtained from high elevation, date the end of MCT thrusting by ca. 15 Ma and record the final stage of exhumation and cooling of MCT termination. All deformation since then is accommodated by the MHT-MBT-MFT fault system. If correct, this implies that the rocks yielding ZHe ages younger than 10 Ma are at least partially or fully reset since the onset of the MBT activity and document the onset of MBT faulting along Himalayan range front in this domain. Previous work showed that across the Pir Panjal Range both thermochronometers have similar age-elevation relationships (Fig. 3). Exhumation rates do not change significantly within the confines of the data resolution. Exhumation is slower, with lower cooling rates in the range interior. Furthermore, this longer transition through the ZHe partial retention and AFT annealing zone may explain some of the observed scatter (Fig. 1C). Deeken et al. (2011) suggested that this exhumation variability within Chamba was best explained as the result of continuous thrusting along a major basal décollement, with a flat beneath the slowly exhuming internal compartments and a steep frontal ramp below the rapidly exhuming frontal range. Their 1D thermal modeling suggests Pliocene–Pleistocene mean erosion and exhumation rates of 0.8–1.9 mm a–1 across the frontal Dhauladhar Range and 0.3–0.9 mm a–1 for the orogen interior Pir Panjal Range over the past ∼15 m.y. Deeken et al. (2011) concluded that the frontal Dhauladhar Range was sufficiently high during this time to build an orographic barrier, focusing climatically enhanced erosional processes and tectonic deformation there. All these results imply a subhorizontal flat northward-dipping MHT below the Pir Panjal Range and no hinterland mid-crustal ramp. However, due to their sparse sampling resolution as well as the strong thermal disturbance caused by rapid lateral and vertical material transport, Deeken et al. (2011) were not able to recognize the main fault accommodating exhumation, ascertain the precise depth of exhumation, or estimate well-constrained exhumation rates.

Pecube-Forward Model Runs

Within the first series of forward model runs, we evaluated four different fault geometry models of the Himalayan wedge for their plausibility (Fig. 4). We compared the model-predicted age patterns with the observed data (Fig. 4). Models III (hinterland ramp in MHT) and IV (duplex) produce a distinct and broad zone of high exhumation rates above the ramp or duplex, resulting in a broad zone of young cooling ages, in contradiction to available age data. For example, models III and IV do not predict mid-Miocene ZHe ages across the Dhauladhar Range and Pir Panjal Range at higher elevations; therefore, these models were not explored further. Model II, with the MHT ramp below the Dhauladhar Range, predicts Quaternary ZHe ages across the entire Dhauladhar Range that are younger than observed, and therefore this model is also dismissed. Only model I, in which the MBT branches off the flat MHT herein, referred to as simple MFT-MBT configuration, predicts cooling age patterns similar to the observations.

For simple MFT-MBT configuration, we explore various geometric and kinematic parameters: (1) variable MBT fault slip rate during period 3 (1–5 mm a–1); (2) variable dip angle of the MBT fault ramp (20°–60°); and (3) variable basal temperature (600–1000 °C) (Fig. S4). We find that 2–4 mm a–1 slip rates along the MBT ramp, 40°–45° dip of the MBT ramp, and basal temperatures of 750–800 °C (Fig. S4c) yield the lowest mismatch between predicted and observed ages using the forward modeling approach. In summary, the forward models suggest that the MBT fault must sole into the MHT to obtain Pliocene reset ZHe ages at the base of southern Dhauladhar Range front. The junction between MBT and MHT must be located ∼30–40 km to the north of the fault trace to reproduce AFT cooling of <3 m.y. invariant with the elevation across the entire Dhauladhar Range (Fig. 1B).

Inverse Modeling

The forward models reveal the primary geometry of the MHT, MBT, and MFT. In order to further constrain the fault slip rates and explore the sensitivity of the model to the assumed thermal parameters, we performed inverse modeling using the Neighborhood algorithm by Sambridge (1999), following the methods presented by Robert et al. (2011), Braun et al. (2012), and Carrapa et al. (2016) to explore the parameter space with lowest misfits. We explored the following four parameters with boundaries chosen to keep a geologically consistent model: (1) the MBT overthrusting rate (1–4 mm a–1); (2) the basal temperature (600–1200 °C); (3) the internal heat production (0–45 °C/m.y.); and (4) the thermal diffusivity (25–60 km2/m.y.). We compare thermochronologic data (αi,dat) and predictions (αi,mod) with a misfit function defined by: 
graphic
where N is the number of data points and σi is the uncertainty of the data. For each inversion, we ran 106–312 models.

The best-fit solutions are insensitive to thermal parameters kept in geologically meaningful ranges (Fig. S6). Predicted cooling ages are controlled mostly by the kinematic field imposed by the MBT slip velocity. For the given MBT geometry (40° dip to a depth of 10 km; 5° dipping MHT), an MBT slip rate of 2.6–3.6 mm a–1 yields the lowest misfit (Fig. 5).

Evaluation of Model Sensitivity

Using best-fit parameters obtained from the numerical inversions (Table 2), we investigate the effect of denudation-induced topographic amplification on the isotherms within the crust. Pecube models allow exploration of the effects of an evolving topography through time through two controls: (1) the topography can be shifted by a vertical offset between two model steps allowing simulation of surface lowering or uplift; and (2) the present-day topography can be multiplied by an amplification factor such that for a factor >1 the relief is amplified, and for a factor <1, the relief is reduced (see Braun, 2003, and Carrapa et al., 2016, for detailed explanations). We explored the effects of variable onset time of MBT activity and of variable evolution of the topography between 10 and 2 Ma using a suite of models designed to simulate different scenarios of faulting and its topographic evolution (Fig. 6).

The first scenario (Fig. 6A) illustrates the effects of variable onset timing on MBT activity and topographic growth during the whole models. We assume a growing topography that begins before the increase of the fault activity. For example, preexisting topography could be explained by previous (early to middle Miocene) shortening across the MCT and therefore could have been established prior to 15 Ma. In all subscenarios, 50% of the present-day topography has been achieved by 10 Ma. When the MBT achieves its full activity, the modern topography is already established.

In the second scenario (Fig. 6B), the growth of topography is induced by the initiation of the MBT faulting. In the different subscenarios, the MBT and topographic growth start together abruptly at different given times after a flat initial topography.

Throughout the third scenario (Fig. 6C), the topographic evolution is decoupled from the onset of faulting. Here we assume that the topography remains constant and equal to the modern topography since the onset of the model.

The last scenario (Fig. 6D) comprises rapid relief growth after the onset of the MBT activity. It accounts for the interior part (Pir Panjal Range) of the model, where the relief growth should decrease through time due to the effective orographic shielding induced by the rising range front.

In summary, all these scenarios reveal that the MBT has to be active for at least ∼8 m.y. to reproduce the observed cooling age pattern, except for the third scenario (Fig. 6C), in which it is achieved after 6 m.y. This value seems to be the response time of the MBT activity to establish the needed thermal field perturbation of the crust to produce the observed age pattern (Deeken et al., 2011). These models show that changes in topography have little impact on the predicted cooling ages. A time-invariant behavior of the predicted ages is achieved after 6–8 m.y., meaning that the landscape is close to exhumational steady state (Willett and Brandon, 2002). The third scenario (Fig. 6C) also predicts ages that are too young. The amplitude of topography needs to have been lower prior to the late Miocene. In contrast, the first and second scenarios yield ZHe ages that are too old. They thus indicate that some topography was probably established prior to the MBT activity and/or fault displacement rates of the MBT were >3 mm/yr.

DISCUSSION

We combined field observations, structural information, thermochronologic data, and thermokinematic modeling to explore ongoing faulting and rapid rock exhumation of the MBT hanging wall along the Dhauladhar Range in the northwest Himalaya. We have conducted sensitivity tests to evaluate minimum onset time of faulting, transferred fault displacement rates into related exhumation rates as well as estimated its contribution to the total shortening across the Himalaya. This setting and these results are strikingly different from most other transects through the MBT further east along strike of the Himalaya.

Field Observations, Horizontal Extent of MBT Fault Ramp, and Their Implications for MBT Activity

In agreement with earlier findings from neighboring sections (Frank et al., 1995), we observed greenschist metamorphic rocks in the hanging wall of the MBT versus unmetamorphosed rocks in the footwall. In the footwall of the MBT, rocks related to the Subathu and Dharamsala Formations are exposed, and are bound to the south by the Palampur thrust (PT, Figs. 1 and 3), which trends parallel to the MBT strike (Fig. 2) (Ranga Rao, 1989; Powers et al., 1998). However, these rocks are internally only slightly deformed, they show no metamorphic overprinting, and they are thrust over the Upper Siwaliks strata. An estimated total fault displacement ranges between 5 and 7 km (the thickness of the Siwalik strata in the region). This is in clear contrast to observations made by Najman et al. (2009), who described the MBT hanging wall at the eastern end of the Kangra reentrant as consisting of unmetamorphosed rocks and illustrating significant lateral variation in the mean fault displacement and exhumation of the MBT hanging wall along strike; their findings were consistent with many other MBT segments (Gansser, 1964; Valdiya, 1980a; Hodges, 2000). In summary, the observed change in metamorphic grade across the MBT fault zone in the Dhauladhar area suggests that within this sector of the northwest Himalaya, the MBT has acted as a primary fault. Furthermore, to a distance of ∼40 km north of the MBT trace, the AFT ages appear to be elevation invariant; this suggests a high vertical rock uplift component and a strong perturbation of its closure isotherms. Such a perturbation of the thermal field likely corresponds to an ∼40 km horizontal length of the MBT fault ramp that merges directly into the MHT north of the Dhauladhar Range (Fig. 1B), forcing the rocks to uplift and cool rapidly with erosion rates high enough to keep up with this. It seems that the topographic wavelength of the Dhauladhar Range exceeds the critical wavelength of the AFT isotherm (Braun, 2002a).

Forward and Inverse Thermokinematic Modeling and Fault Geometries

We tested several fault geometries of the northwestern Himalayan wedge using forward thermal models. Only scenario I (Fig. 4), where an imbricated MBT fault ramp roots into the MHT, could explain the observed cooling pattern derived from the low-temperature thermochronology data (geometry I, Fig. 4). Neither a hinterland MHT ramp geometry (scenario III, Fig. 4) nor a duplex beneath the Pir Panjal Range was sufficient to reproduce the observed cooling age pattern (scenario IV, Fig. 4). The exhumation rates obtained from such geometries are too high and they yield model ages that are too young. Therefore, we interpret the Dhauladhar Range as the result of continuous slip along the MBT ramp.

We systematically tested a range of possible values of dip angles (20°–60°) of the MBT fault ramp using forward models (S3 and S4 in the Data Repository). We found that dips of the MBT ramp between 35° and 45° yield the lowest misfit. These dips are in agreement with our own field measurements and Adlakha et al. (2013). In contrast, DeCelles et al. (2001) inferred a ∼10° dip of the MBT in the central Nepal segment.

Our detailed inverse models suggest uniform slip rates of ∼2.6–3.6 mm a–1 along the MBT (Fig. 5; Table 2) over at least the past 6 m.y. Together with our assumed structural and thermochronologic constraints on the estimated duration of faulting (∼10 m.y.), the minimum total fault displacement on the MBT is in the range of 26–36 km. At the position where the MBT soles into the basal décollement, the thermal gradient obtained from our numerical inversions would result in a temperature slightly >300 °C. Assuming a thermal gradient of ∼30 °C/km, this is a depth of ∼10 km and therefore agrees with the depth of the MHT (sole thrust) and the structural architecture defined in our numerical models. Other geometries would exhume material from greater depth and result in exposure of higher metamorphic facies at the surface, as observed in the field and described herein.

The observed garnet growth in the MCT hanging wall (Fig. 2), however, does not fit with our modeling. Several possibilities could explain the existence of the garnets: (1) they could have grown at higher pressure-temperature conditions at deeper depths and were transported first southward along the décollement before they were uplifted across the MBT ramp; or (2) the temperature conditions could have been higher due to increased heat advection along the MHT shear zones (e.g., due to fluid circulation); (3) the fact that we found these garnets so close to the MCT fault zone could also hint that Higher Himalayan material has been transported by the MCT from a far interior source during the early to mid-Miocene. In summary, the inverse models are in good agreement with intense and long-lasting exhumation above a steep frontal MBT fault ramp (Fig. 1).

We tested the sensitivity of the models to different MBT initiation ages (see Fig. 6). These models show that a minimum duration of activity for the MBT of at least 6 m.y. is necessary to reproduce the observed cooling age pattern. However, our results are not strongly sensitive to changes in topography or initial topography. At least 6–8 m.y. are needed to establish a time-invariant pattern for AFT closure isotherms. Such patterns are established when exhumation is in a near thermal steady state (Willett and Brandon, 2002). Our results document that local exhumation history is strongly interconnected with the underlying fault geometry. The most realistic exhumation behavior occurs when the Dhauladhar Range reaches a threshold elevation (>3 km), where it forms an orographic barrier. Rock uplift is then balanced by rapid erosion across the entire Dhauladhar Range (Fig. 6D; Deeken et al., 2011). Overall, our model results are in good agreement with sedimentary constraints suggesting a time of onset of MBT activity by 10–8 Ma (Meigs et al., 1995; Brozovic and Burbank, 2000). Assuming that the local climate pattern has been approximately consistent since the Dhauladhar Range was established, this range formed the primary orographic barrier soon after ca. 10 Ma, as the present-day precipitation pattern documents (Bookhagen and Burbank, 2006). A possible climatic feedback between rapid erosion and denudation-induced unloading and rock uplift (Wobus et al., 2003; Thiede et al., 2004) would maintain the steep southern flanks of the Dhauladhar Range and explain the rapid exhumation rates and localization of the deformation over the past 10 m.y. Such a time frame is unusual compared to other MBT segments along the entire Himalaya, but agrees with findings by Deeken et al. (2011). Therefore, the Dhauladhar Range may have achieved exhumational steady state for at least the last couple of million years. No seismicity or geomorphic evidence, such as fault scarps or offset fluvial terraces, related to any present-day MBT activity have been recognized. This is probably because vigorous hillslope and fluvial erosion processes along a steep mountain front could rapidly erase any exposure of transient sedimentary deposits or small tectonic landforms. Nevertheless, our thermokinematic modeling indicates that the MBT is potentially still an active fault. Based on the mean fault displacement rates (2.6–3.6 mm a–1) and best-fit MBT dip, this implies that the MBT within this domain has probably accommodated ∼1.7–2.4 mm a–1 of crustal shorting over the past 10 m.y. Assuming that present-day geodetic shorting rates (∼14 ± 2 mm a–1) across the northwest Himalaya have been sustained over geologic time scales, this would imply that here the MBT accommodated ∼15% of the total Himalayan convergence from its onset. This has not been recognized by previous workers, neither by various geodetic studies due to insufficient data resolution nor by field studies, which favor an accommodation of the total Himalayan crustal shortening since 2 Ma within the Sub-Himalaya (Powers et al., 1998; Thakur et al., 2014).

Overall, our data support the view that the MHT décollement dips uniformly ∼3°–5° to the northeast and indicate an absence of any hinterland mid-crustal ramp along the MHT within this segment, in contrast to other Himalayan segments further east (Valdiya, 1980b; Srivastava and Mitra, 1994; Webb, 2013). The obtained cooling pattern signifies a clear relationship between morphology and underlying fault geometry and emphasizes the outstanding role of the Dhauladhar Range compared to the Himalayan front further east.

CONCLUSIONS

Our new results along the southern flanks of the Dhauladhar Range strongly suggest that within this sector of the northwest Himalaya, the MBT is a primary fault and possibly still active. This is supported by the observed change in metamorphic grade across the MBT fault zone. Combining new and previously published ZHe and AFT data along with Pecube 3D-thermokinematic modeling of the Dhauladhar Range provides key information about the exhumation history of the range front. Our key findings are the following.

From our low-temperature thermochronologic data, we infer that the MBT, in contrast to other Himalayan segments, forms a steep fault ramp that soles into the gently dipping MHT décollement and has exhumed rock from depths of 8–10 km since the late Miocene.

Our modeling results reveal that the ZHe and AFT cooling ages can only be explained with fault slip at mean rates of ∼2.6–3.6 mm a–1 since the MBT activation and along a 30–40-km-long MBT fault ramp that directly soles into the MHT. Since the late Miocene, the MBT has accommodated ∼1.7–2.4 mm a–1 horizontal shortening or nearly 15% of the total Himalayan crustal shortening across this sector.

These results imply strong along-strike structural variations in the fault geometry of the Himalayan orogenic wedge. The Chamba segment exhibits strikingly dissimilar structural pattern compared to the central Himalaya. For example, the cooling age pattern of the Chamba transect is explained by faulting along flat or gently northeast dipping MHT décollement. It excludes a hinterland mid-crustal ramp or the duplexing of the Lesser Himalayan sequence in this sector of the range.

ACKNOWLEDGMENTS

We appreciate the collaboration and logistical support provided by V. Jain and his group (Indian Institute of Technology Gandhinagar, Ahmedabad). This study was funded by a Deutsche Forschungsgemeinschaft (DFG) grant to Thiede (TH 1371/5–1) and a DAAD-DST (Deutscher Akademischer Austauschdienst–India Department of Science and Technology; public-private partnership PPP-India #57035520) grant to Thiede and V. Jain. Faruhn and Dey received support from DFG-GRK (Graduiertenkollegs) grant 1364 to M. Strecker (STR 373/19–2). We thank T. Ehlers for collaboration and ZHe-dating support at the University of Tubingen, P. van der Beek for hosting Thiede at ISTerre (Institut des Sciences de la Terre), Université Grenoble Alpes, and R. Arrowsmith for comments on an early version of this manuscript. We also thank Birgit Fabian for support with figure drafting. The comments of two anonymous reviewers are greatly appreciated and helped to improve the manuscript significantly.

1GSA Data Repository Item 2017241, which contains Supplemental figures S1–S6 and one supplemental table providing the analytic details about ZHe dating and about the Pecube modelling and data applied, is available at http://www.geosociety.org/datarepository/2017, or on request from editing@geosociety.org.