Tectonic displacement of drainage divides and the consequent deformation of river networks during crustal shortening have been proposed for a number of mountain ranges, but never tested. In order to preserve crustal strain in surface topography, surface displacements across thrust faults must be retained without being recovered by consequent erosion. Quantification of these competing processes and the implications for catchment topography have not previously been demonstrated. Here, we use structural mapping combined with dating of terrace sediments to measure Quaternary shortening across the Indus River valley in Ladakh, NW Himalaya. We demonstrate ∼0.21 m k.y.–1 of horizontal displacement since ca. 45 ka on the Stok thrust in Ladakh, which defines the southwestern margin of the Indus Valley catchment and is the major back thrust to the Tethyan Himalaya in this region. We use normalized river channel gradients of the tributaries that drain into the Indus River to show that the lateral continuation of the Stok thrust was active for at least 70 km along strike. Shortening rates combined with fault geometries yield vertical displacement rates that are compared to time-equivalent erosion rates in the hanging wall derived from published detrital 10Be analyses. The results demonstrate that vertical displacement rates across the Stok thrust were approximately twice that of the time-equivalent erosion rates, implying a net horizontal displacement of the surface topography, and hence narrowing of the Indus Valley at ∼0.1 m k.y.–1. A fill terrace records debris-flow emplacement linked to thrust activity, resulting in damming of the valley and extensive lake development. Conglomerates beneath some of the modern alluvial fans indicate a northeastward shift of the Indus River channel since ca. 45 ka to its present course against the opposite side of the valley from the Stok thrust. The integration of structural, topographic, erosional, and sedimentological data provides the first demonstration of the tectonic convergence of drainage divides in a mountain range and yields a model of the surface processes involved.
The topography of active mountain ranges records surface uplift in response to crustal thickening countered by erosion (e.g., Dahlen, 1990). The horizontal velocities that drive crustal thickening are commonly an order of magnitude higher than the vertical velocities, and so it is expected that this should be recorded by the topography (Pazzaglia and Brandon, 2001; Willett et al., 2001; Miller and Slingerland, 2006). Model experiments have indicated that the broad asymmetry of many small mountain ranges, such as the Southern Alps of New Zealand, the Pyrenees, and ranges in Taiwan, may be explained by the horizontal translation of deforming rock from the side of the range dominated by accretion toward the opposing side (Willett et al., 2001; Sinclair et al., 2005; Herman and Braun, 2006).
It is reasonable to suggest that such large-scale forcing of topography must also play a role in determining the geometry of river catchments and their channel courses. At the largest scale, it is proposed that the extraordinarily elongate forms of the rivers draining eastern Tibet (Salween, Mekong, and Yangtze) represent highly strained forms of previously more regularly shaped catchments in response to distributed crustal shortening and rotation around the eastern corner of the Indian indentor (Hallet and Molnar, 2001). Similarly, the river catchments of the Southern Alps of New Zealand are understood to have been deformed to their present shape during oblique convergence (Koons, 1995; Castelltort et al., 2012). Tectonically induced changes in catchment shape may be further modified by river capture and progressive migration of drainage divides in response to factors such as variability in rock strength (Bishop, 1995), changing river base levels (Mudd and Furbish, 2005), and ridge-top glaciation (Dortch et al., 2011a). The competition between tectonic deformation of river catchments and the response of the rivers is highlighted across the Himalaya, where all of the big rivers are characterized by steepened reaches and more localized knick zones as they respond to variable rock uplift fields (Seeber and Gornitz, 1983; Wobus et al., 2006a). The smaller river catchments near the foothills of the Himalaya exhibit variable catchment geometries in response to lateral advection over thrust ramps (Champel et al., 2002; Miller et al., 2007). Large-scale catchment deformation has broad implications for the topographic form of active mountain ranges and the distribution of erosion and transported sediment to surrounding sedimentary basins. Any modification of catchment shape also has implications for the scaling of upstream catchment area with channel length and hence the long profile of rivers (Whipple and Tucker, 1999; Willett et al., 2014).
Fluvial erosion into bedrock can be approximated by a power-law relationship between channel slope and river discharge (Howard et al., 1994; Whipple and Tucker, 1999). In this stream power model, the fault offset generates an oversteepened channel reach (knickpoint or knick zone) that migrates upstream as a kinematic wave. Additionally, the model predicts that sustained differential rock uplift across a fault will generate increased channel steepness (for a given upstream area) on the upthrown block. Analysis of channel steepness has been used to assess fault activity in mountain ranges (e.g., Hodges et al., 2004; Kirby and Whipple, 2012), with relative rock uplift in the hanging wall of a thrust fault leading to increased stream power generated by channel steepening.
Little is known of the interaction between thrust shortening and the consequent deformation of catchment shape, as opposed to the offset of individual channels by faults. As yet, there has been no demonstration of the horizontal convergence of drainage divides in response to shortening on a thrust fault that bisects a catchment. The challenge that this sets is the requirement to quantify both the shortening and the time-equivalent erosional response.
The objective of this study is to test whether rates of horizontal displacement across a thrust fault are capable of driving the horizontal convergence of opposing drainage divides when moderated by the erosional response to fault displacement. We examine the Indus River valley in Ladakh, NW India, which is one of the largest longitudinal river catchments of the Himalaya, with an average width of around 35 km and a length of ∼200 km parallel to the mountain range in this region. The aim is first to test for the presence of active shortening across the Indus Valley, which has never been demonstrated. This is regionally significant because the valley follows the line of the main back thrust in the region, carrying Tethyan Himalaya units northeastward toward the Gangdese batholith (van Haver, 1984; Searle et al., 1990). Large portions of the Indus and Tsangpo Rivers further east in the Himalaya also follow this structural feature. Thrust displacement rates were measured using mapped and dated alluvial and lacustrine terraces, and by documenting displacement of these terraces across faults. Second, having presented evidence for Quaternary deformation, we compare the vertical component of rock displacement in the hanging wall of the main back thrust relative to the magnitude of erosion at similar time scales (Fig. 1); this is this ratio that determines the signal of topographic change across the valley. Erosion rates are presented using published low-temperature thermochronology (Kirstein et al., 2006, 2009) and detrital cosmogenic nuclides (Dortch et al., 2011a; Munack et al., 2014; Dietsch et al., 2014). In addition, the distribution of changing erosion rates in response to thrust displacement is inferred regionally through an analysis of river channel steepness of catchments that drain into the Indus Valley. Third, sedimentological evidence for valley damming in response to fault movement, and for the migration of the main Indus channel, is presented. The integration of structural, topographic, erosional, and sedimentological data provides the first quantification of the tectonic convergence of drainage divides in a mountain range and yields a model of the surface processes involved.
The Indus River of Ladakh flows northwestward (Fig. 2) between the highly deformed Cretaceous to Miocene sediments of the Indus molasse, which are thrust northeastward against the relatively undeformed Cretaceous and Paleogene Ladakh batholith complex (Figs. 2 and 3). The Indus molasse records sedimentation in a forearc basin that evolved into an intramontane basin following continental collision (Garzanti and Van Haver, 1988; Searle et al., 1990; Sinclair and Jaffey, 2001). The Ladakh batholith forms part of the Gangdese batholith complex at the boundary between the northern mountains of the Himalaya and the Tibetan Plateau. It represents the magmatic arc prior to continental collision and consists of a succession of granodioritic rocks overlain by a volcanic succession that forms the southern wall of the Shyok Valley to the north (Weinberg and Dunlap, 2000).
The Indus molasse of the Stok Range is intensely deformed, with fold-and-thrust structures verging to the northeast and southwest. At the boundary with the Ladakh batholith, the Cretaceous succession locally onlaps the margin of the batholith (Van Haver, 1984), but the main topographic boundary is defined by a thrust fault that carries steeply tilted Miocene molasse successions in its hanging wall over Quaternary alluvial-fan deposits; we term this the Stok thrust (Fig. 3), which laterally correlates to the Great Counter thrust further east (Murphy and Yin, 2003). The bulk of deformation of the Indus molasse has occurred since deposition of the youngest sediments around 20 Ma (Sinclair and Jaffey, 2001). The extent to which deformation has continued since this time has not been documented.
The Ladakh batholith contains crystallization ages ranging from ca. 103 to 47 Ma (Honegger et al., 1982; Weinberg and Dunlap, 2000), and it is overlain by a volcanic succession along its northern margin that is tilted steeply northeastward (Weinberg et al., 2000). This rotation is thought to have occurred in the hanging wall of a thrust fault that dips northeastward under the batholith, and which was active during early Miocene times (Kirstein et al., 2006); this structure is comparable to the Gangdese thrust near Lhasa (Yin et al., 1994). Thermochronological analyses using apatite and zircon U-Th/He dating and apatite fission-track dating indicate rapid cooling of ∼25 °C/m.y. around 22 Ma, followed by a deceleration to rates <3.5 #xb0;C/m.y. since then (Kirstein et al., 2006).
Detrital cosmogenic 10Be analysis across the Ladakh batholith indicates erosion rates of ∼0.04–0.09 m k.y.–1 for the main tributaries on the northeastern side and 0.02–0.05 m k.y.–1 on the southwestern side (Dortch et al., 2011a; Munack et al., 2014). Smaller, side tributaries on the southwestern side of the batholith record rates as low as 0.008 m k.y.–1 (Dietsch et al., 2014); these represent the slowest rates recorded from the Himalaya. These measurements average over tens of thousands of years and record an asymmetry in erosion rates associated with greater degrees of glaciation on the northern side of the Ladakh batholith driving glacial headwall erosion and migration of the drainage divide toward the southeast over this time period (Jamieson et al., 2004; Dortch et al., 2011a). Small (∼1-km-long) glaciers are still present at the drainage divide around 5500 m elevation, with significant glacial erosion having occurred down to ∼4700 m on the southwestern side of the batholith (Hobley et al., 2010). Dating of boulders on moraines in the Ladakh region has demonstrated multiple glaciations in this region (Owen et al., 2006; Owen and Dortch, 2014), with the oldest significant glaciation being ca. 300 ka (Dortch et al., 2013). On the southwestern margin of the Indus Valley, erosion rates from the Indus molasse successions of the Stok Range are faster than on the batholith, with 10Be concentrations implying millennial erosion rates of 0.07–0.09 m k.y.–1 (Munack et al., 2014).
In the Leh region of the valley, the northwesterly flowing Indus River is bounded by large alluvial fans draining the Indus molasse from the southwest. These fans appear to force the present river channel to bank up against the interfluve ridges of the batholith to the northeast (Fig. 3). A terrace containing evidence of lake sedimentation forms the distal margin of these alluvial fans (Fig. 4), and other terraces in the valley testify to a history of damming of the Indus River (Bürgisser et al., 1982; Fort, 1983; Phartiyal et al., 2005; Blöthe et al., 2014). The presence of broad regions of alluvium in the lower reaches of the tributaries draining the batholith (geomorphic domain 3 of Hobley et al., 2010, 2011) encouraged Jamieson et al. (2004) to suggest that an asymmetry in deformation and erosion across the Indus Valley has resulted in northeastward translation of the valley over the batholith. However, evidence for ongoing structural deformation and relative displacement of the Indus molasse has not been recorded (Dortch et al., 2011a) and is therefore a key focus of this study. As the valley is traced northwestward from the village of Phey, so the river’s course cuts a large gorge into the deformed molasse, and the long profile exhibits a broad steepening downstream of the alluviated reach in the Leh Valley (Jamieson et al., 2004).
EVIDENCE FOR QUATERNARY SHORTENING
Fan Terrace Data
Geomorphic fill terraces usually record abandoned floodplain surfaces that parallel the modern river channel, and they can usually be correlated across the landscape, and so can be used to assess evidence of ongoing deformation since formation (e.g., Lavé and Avouac, 2001; Pazzaglia and Brandon, 2001; Wegmann and Pazzaglia, 2009). The terraces in the Leh region of the Indus Valley represent the abrupt downslope termination of alluvial-fan surfaces into a 20–80 m succession of bedded sandstones and laminated siltstones that record floodplain and shoreline settings around the edge of ancient lakes; this sedimentological transition is associated with a geomorphic break recording the approximate coastline of the paleolake. In order to distinguish these features from classic fill terraces (e.g., Wegmann and Pazzaglia, 2009), we refer to these as “fan terraces.” One of the best documented sections through a fan terrace succession is in the Spituk region near Leh, where radiocarbon dates yield ages from ca. 51 to 31 ka (Phartiyal et al., 2005). Several terrace successions also contain extensive soft sediment deformation that has been interpreted as a record of seismicity throughout the region (Phartiyal and Sharma, 2009). We mapped two terrace fill successions around the northwestern part of the Leh Valley that could be correlated across the two sides of the valley (Fig. 4). Field mapping of terrace successions using a laser range finder was supported by Google Earth satellite imagery and the 1 arc second Shuttle Radar TopographyMission digital elevation model (DEM); 1 arc second equates to approximately 30 m horizontal resolution. The top surface of the higher fan terrace (T1) is at an average elevation of around 3250 m and represents the dissected remnant of an alluvial fan, with lacustrine sediments at a downslope break in topography (Figs. 4 and 5). A lower fan terrace succession is capped by a surface (T2) at around 3200 m elevation and is evident throughout the region. This level forms the break of slope between alluvial fans that drain the Stok Range and the modern Indus River floodplain in the Leh Valley (Fig. 4).
The sedimentology of the T1 infill is best exposed around Spituk (Fig. 4), where at least 50 m of silts, sands, and gravels are present (Fig. DR11), recording marginal lake sedimentation (Bürgisser et al., 1982; Phartiyal et al., 2005). The lower portions of the section are dominated by coarse-grained, fining-upward event beds delivered from marginal deltaic feeder systems. This thick succession underlying the T1 surface can be traced at the same elevation downstream for at least 10 km (Fig. 4B). The lower T2 infill is exposed in the cliffs on the southwest side of the valley opposite Spituk. This succession is ∼20 m thick and dominated by poorly bedded coarse gravels and breccias typical of alluvial-fan sedimentation. Approximately 2–4 m below the fan surface, there is a succession of well-bedded, fine to medium sands with some planar lamination, and some evidence of rootlets, grass blades, shells, and other organic material. There is also a 40 cm unit of finely laminated siltstones, similar to the lacustrine deposits of the T1 fill (Fig. DR2 [see footnote 1]). This interval is interpreted as an episode of lacustrine and marginal floodplain sedimentation that defined the base level for the alluvial fans that drain the Stok Range (Fig. 4A). In contrast to the T1 fill succession, downstream tracing of the T2 terrace fill demonstrates a reduction in elevation that is parallel, but ∼25 m above the modern Indus River.
As the Indus River continues downstream to the northwest, so it changes course from flowing at the boundary between the Indus molasse and the Ladakh batholith to flowing within, and along the strike of the Indus molasses, where it forms a steep gorge (Figs. 4 and 5A). On either side of this gorge, the two terraces fills are clearly visible, with the T1 fill characterized by light, cream-colored lake sediments, and T2 having a more pink tone where the sediment forms a bench in the gorge. Near to the turning for the Markha Valley, the T1 fan terrace is deformed by thrusting, folding, and extensive irregular soft sediment deformation (Fig. 6). At the southwestern extent of the terrace, it is overthrust by the Indus molasse on a fault dipping 37° to the southwest. Thinly bedded alluvium is folded into a broad syncline in the footwall of the fault with a wavelength of ∼200 m (unit 1; Figs. 6 and 7). Within this lower succession, there are meter-scale thrust faults and folds that are draped by overlying beds and hence are syndepositional. An unconformity divides this folded succession into two, recording a phase of erosion and renewed sedimentation prior to the final phase of folding. These folded alluvial sediments are truncated by a structureless breccia with meter-scale blocks of the Indus molasse that is interpreted as a surficial debris-flow deposit that ranges from 2 to 5 m thick (unit 2; Figs. 6 and 7). This debris flow is draped by finely laminated pale siltstones that are interpreted as lake deposits (unit 3; Figs. 6 and 7). These siltstones are capped by gravels of the abandoned T1 alluvial-fan surface (unit 4; Figs. 6 and 7); this surface has since been dissected by a dense network of modern river channels (Fig. 5B). In comparison, the lower T2 fill is undeformed.
These exposures are interpreted as a record of syndepositional thrust faulting that caused progressive deformation of Indus Valley alluvium. Coseismic slope instability resulted in extensive landslide material, which was subsequently reworked during a storm to deposit a debris flow that dammed the valley, leading to lake formation. Folding and intraformational unconformities in the footwall of the thrust indicate that this records fault propagation folding with associated growth strata (e.g., Suppe et al., 1992). Soft sediment deformation has been recognized elsewhere in this T1 terrace fill, as well as a fault offset between the batholith granites and lake sediments near Spituk (Phartiyal and Sharma, 2009).
The exposures in the region of Spituk, near Leh (Fig. 4), were dated using four radiocarbon ages that range from ca. 50.8 ± 5 ka at the base to ca. 31.0 ± 0.7 ka near the top (Fig. 3; Phartiyal et al., 2005). Given the significance of thrust shortening of the T1 terrace, we chose to date the deformed T1 terrace sediments near the Markha Valley using optically stimulated luminescence (OSL) on quartz, and to compare this against the range of radiocarbon ages at Spituk, and against new ages for the other terraces.
We collected 20 samples of medium- and fine-grained sand and silt layers for OSL dating of quartz grains, and most samples were derived from units that were interbedded with coarse-grained or conglomeratic deposits of fluvial and alluvial-fan origin (see supplementary material for full description [see footnote 1]). Other deposits that were sampled record lacustrine environments, and reworked horizons overlying mass-flow deposits. Samples were collected in copper tubes (2.5 cm diameter, 12 cm long) that were tapped into the target deposits parallel to the stratigraphic orientation. The tubes were sealed with black tape to avoid light penetration and to minimize any moisture loss within the tubes. At least 2 cm of sediment from both ends of each tube were used for dosimetry measurements, and the remaining material was used for dating. Analysis of luminescence behavior, dose rate estimation, and age calculations were conducted at University of St. Andrews using the protocol outlined in King et al. (2013). The analytical details and results (with tables and figures) are presented in the supplementary material (see footnote 1). Only 17 of the 20 samples were dated, and two ages are based on a low number of aliquots (Zansk2011-1 and Nimmu2011-1).
T1 terrace fill. The four samples from the deformed T1 terrace succession near the Markha junction generated ages, in ascending stratigraphic order, of 35.6 ± 2.7 ka, 73.0 ± 0.7 ka, 40.0 ± 5.2 ka, and 77.2 ± 11.7 ka (Figs. 6 and 7); given the observed stratal sequence, these cannot all represent true depositional ages. Having confidently correlated the T1 succession from Spituk to the Markha junction (Fig. 4B), we would expect the ages to fall within the time interval of 50.8 ± 5.4–31.0 ± 0.7 ka based on the radiocarbon ages at Spituk (Phartiyal et al., 2005). In order to be confident of the OSL correlation to the radiocarbon ages, we also ran a sample from the top of the Spituk T1 succession and obtained an age of 27.5 ± 3.0 ka, which is within error of the radiocarbon age. Consequently, our interpretation of the ages at the Markha junction locality is that the two ages that are significantly older than the radiocarbon age bracket record age overestimation. Inheriting older ages is common in fluvial systems where sediment grains were not fully exposed during transport and deposition, meaning that their luminescence “clocks” had not been reset (incomplete bleaching; Wallinga, 2002). This is particularly common where coarse sands were deposited by short-lived, turbid flows and mass flows, which are typical in alluvial-fan settings.
T2 terrace fill. The T2 terrace fill was sampled on the opposite side of the valley from Spituk at the margins of the large alluvial fans that dip gently northward into the Indus Valley (Fig. 4). The samples were taken ∼20 m above the modern floodplain and consisted of sands and gravels with finer-grained intervals (Fig. DR2 [see footnote 1]). The three samples (Dung2011-01, Dung2011-02, and Dung2011-03) yielded ages of 22.0 ± 1.3 ka, 19.1 ± 0.7 ka, and 11.7 ± 0.7 ka (Table DR3 [see footnote 1]). Similar ages ranging between 22.0 ± 1.3 ka and 8.8 ± 0.8 ka from terrace levels downstream near Nimu and Basgo suggest that this was a period of widespread sediment aggradation throughout this part of the Indus Valley.
Horizontal Displacement Rates
Based on the stratigraphic location of the T1 Markha junction samples (Fig. 6B), deformation of this succession must have started during accumulation of the alluvial deposits of unit 1 with ages of 35.6 ± 2.7 ka and 40.0 ± 5.2 ka; in order to convey conservative estimates of shortening rates, we use the oldest possible age for unit 1 of 45.2 ka. The horizontal shortening at the Markha junction locality was measured from three components of the deformation (Fig. 6C): (1) the minimum value of 2.3 m of horizontal displacement of the Indus molasse at the thrust (the molasse is not observed in the footwall of this thrust, so the value could be much greater); (2) a line-length restoration of the folding in the lowermost unit 1, beneath the unconformity, amounting to 4.5–5.5 m of shortening; and (3) a line-length restoration of the folding of the upper part of unit 1, totaling 2.0–2.5 m, which was then truncated by the debris flow (Fig. 6C). The errors in these line-length restorations are based on uncertainties in projected fold limbs. By combining the midpoint values of these three shortening measurements, we obtain a total shortening of at least 9.5m. Based on the total horizontal displacement (folding and faulting) since deposition of unit 1 of 9.5 m, and the oldest age for the deformed alluvium of 45.2 ka, we estimate a mean shortening rate from that time to the present of at least 0.21 m k.y.–1.
Topographic Expression of Shortening
Whether the activity on the Stok thrust was localized or regional is significant in the context of its impact on orogenic topography. Therefore, we used fluvial topography to test the lateral extent of thrust activity in the Indus molasse (e.g., Kirby and Whipple, 2012).
It has been recognized for over a century that erosion rates in bedrock channels should increase with increasing channel gradient and water discharge (e.g., Gilbert, 1877). If other factors are equal, for example, rock hardness or local uplift rates, channel gradients should decrease as discharge (or its proxy, drainage area) increases, and so any topographic analysis that uses channel gradients as a proxy for erosion rates must take into account drainage area. Several authors have used a scaling relationship to explore changes in erosion rates along bedrock channels: S = ksA–θ, where S is the topographic slope, ks is a steepness index regressed from slope and area data, A is the drainage area, and θ describes the rate of change of slope or concavity of the long river profile, (Wobus et al., 2006b). If θ is set to a fixed value, the steepness index ks becomes a normalized steepness index, ksn, and this index has been applied to a number of regions of active tectonics; importantly, it can identify differential rock uplift fields that are bordered by faults that have not been historically active, and so aid seismic hazard awareness (Kirby and Whipple, 2012). However, the selection of θ and identification of reaches with statistically different values of ksn can be difficult with noisy slope and area data.
Our topographic analysis of river long profiles normalizes for drainage area by integrating drainage area over flow distance. This method, first suggested by Royden et al. (2000), produces a transformed coordinate, χ (chi), which has dimensions of length (Perron and Royden, 2013). The elevation of the channel can then be plotted against the χ coordinate, and the gradient of the transformed profile in χ-elevation space provides a steepness indicator that can be used to compare channel segments with different drainage areas.
The transformed coordinate is calculated with
where x [length] is the flow distance from the outlet, with dimensions henceforth denoted as length [L] and time [T] in square brackets, xb [L] is the flow distance at the outlet, A [L2] is the drainage area, A0 [L2] is a reference drainage area introduced to ensure the integrand is dimensionless, and m and n are empirical constants, and where –m/n = θ.
where E [L T–1] is the erosion rate, S [dimensionless] is the slope, and K is an erodibility coefficient with dimensions that depend on the exponent m. Royden and Perron (2013) demonstrated that in landscapes where channel incision could be described by Equation 2, changes in erosion rates at the base of channels would result in upstream-migrating “patches” or segments of constant slope in χ-elevation space, given constant bedrock erodibility and local uplift rates. These segments can be described by:
In channels with fluvial incision that can be described with the stream power law (Eq. 2), the gradient in χ-elevation plots is closely related to the normalized channel steepness (ksn), which is simply the steepness index (ks) that has been calculated using a reference concavity, θ (e.g., Hodges et al., 2004; Kirby and Whipple, 2012):
Other models have been proposed for channel incision, including those that incorporate the role of sediment supply (Sklar and Dietrich, 1998) and erosion thresholds (e.g., Snyder et al., 2003). However, even if the stream power incision model is an imperfect description of channel incision (cf., Lague, 2014), Gasparini and Brandon (2011) demonstrated that Equation 2 works as an approximation of the proposed incision model. At a minimum, both Mχ and ksn can still be calculated and allow a qualitative comparison of the steepness of channel segments relative to their upstream area from different parts of the channel network. Both chi-analysis and the normalized steepness index (ksn) have been found to correlate well with erosion rates in the Yamuna River, which is a basin to the south of Ladakh (Scherler et al., 2014), and many other studies have shown a correlation between steepness index and erosion rate (Ouimet et al., 2009; DiBiase et al., 2010; Miller et al., 2013; Kirby and Whipple, 2012; Harel et al., 2016).
We use a method developed by Mudd et al. (2014) to determine the most likely locations of channel segments, defined as sections of the channel that have distinct slopes in chi-elevation space. Our method is aimed at using channel profiles to infer differences in erosion rates and thus has the same objectives as previous studies that used (ksn). The chi method differs from previous studies based on slope-area information in that it uses elevation rather than slope, which is less noisy (e.g., Perron and Royden, 2013).
The Mudd et al. (2014) method tests all possible contiguous segments in a channel network and selects the most likely segment transitions using the Akaike information criterion (AIC; Akaike, 1981), which is a statistical technique that rewards goodness-of-fit while at the same time penalizing overfitting. Mudd et al. (2014) used both field examples and numerical models to show the method could distinguish channel segments of varying erosion rates via detection of varying Mχ values; their results followed the analytical work of Royden and Perron (2013), demonstrating that the chi method could distinguish varying erosion rates in transient landscapes. Changes in Mχ may be due to factors other than changing erosion rates; for example, changes in channel erodibility could force changes in Mχ. The Mudd et al. (2014) method is agnostic with regards to the cause of changing Mχ values, it simply finds segments with different Mχ values that may be differentiated statistically.
Because channel steepness using Mχ or ksn is viewed through the lens of a detachment-limited incision model (i.e., Eq. 2), channels compared with these methods should be detachment and not transport limited. Channels to the south of the Indus, upstream of the fans, are bedrock, but channels to the north of the Indus contain coarse sediment (Hobley et al., 2010). However, Hobley et al. (2011) demonstrated through both field measurements and numerical modeling that the evolution of these channels could be best described with a detachment-limited model of channel incision, and should therefore be amenable to Mχ or ksn analysis.
To calculate both segments and Mχ values, the transformation of Equation 1 requires values for both A0 and m/n to be selected. The reference drainage area simply scales χ, so it changes the absolute of Mχ but not relative values. The m/n ratio is, on the other hand, determined statistically, following the method of Devrani et al. (2015), in which target basins are selected (in our case, 14 basins), and in each basin, 250 sensitivity analyses are run in each of the 14 basins to determine the range of m/n ratios in each basin and to determine a regional m/n value to be used in calculating Mχ values. Of the 14 basins, seven were located in the Ladakh batholith, and seven in the molasse, and in each, we ran the collinearity algorithm of Mudd et al. (2014). This algorithm both selects the most likely segments in the channel network and also seeks to maximize the collinearity of the main stem and tributaries, which is expected if the most plausible m/n is selected (Perron and Royden, 2013). The Mudd et al. (2014) algorithm takes several parameters: the minimum length of a segment, the maximum number of nodes in a subsection of profile to be analyzed, the uncertainty of the elevation, and a parameter that determines how frequently sections are sampled by a Monte Carlo sampling routine (each m/n analysis is the result of several hundred resampling runs of the channel network). The m/n ratio can change depending on these parameters, so we ran 250 iterations across parameter values and selected the most likely m/n ratio. Thus, the m/n ratio was selected based on 250 × 250 × 16 implementations of the Mudd et al. (2014) algorithm.
In the Ladakh region, it has been previously noted that river concavity (i.e., m/n) varies in relation to the degree of upstream glaciation (Hobley et al., 2010), which suggests that local channel slopes are not a simple function of rock uplift or lithology. In order to avoid the influence of glaciation, we selected those catchments where moraines, valley widening, and channel slope reduction due to glacial erosion were absent—these being the characteristics of the upper glaciated domain of Hobley et al. (2010). Based on 14 of these smaller, nonglaciated catchments (Fig. 8B), we derived a range of concavity values with a mean of 0.4 for both the Indus molasse and the Ladakh batholith; i.e., there was no significant difference between them. Once we determined the regional m/n ratio, we then applied this to all the river networks in the region to map Mχ values of channels draining into the Indus from both the north and south.
The channel steepness for all rivers across the region demonstrated a high degree of variability (Fig. 8A), particularly within the larger tributaries that drain from the glaciated drainage divide of the Ladakh and Stok Ranges. These variable Mχ values link directly to the three geomorphic domains associated with glacial erosion, incision into glacial moraines, and alluvial-fan growth identified by Hobley et al. (2010). Therefore, these larger catchments were not used for the evaluation of variable erosion rates across the region; we speculate that the variation may be linked to the sediment flux–dependent channel incision processes documented in a number of these valleys by Hobley et al. (2011). However, the smaller unglaciated catchments, which range from 4 to 18 km in length, provide Mχ values that can be compared throughout the region (Fig. 8B).
We compared Mχ values for opposing catchments on either side of the Indus River valley (Fig. 9C), which, due to the proximity of their outlets, have the same local base level (i.e., the Indus River). Mχ values are consistently higher, and more variable, on the southwestern margin of the valley on the Indus molasse compared to the opposing tributaries that drain the batholith. Within the batholith, the relatively constant Mχ values suggest that there is little spatial variation in local uplift rates, channel erodibility, or erosion driven by base-level changes; it is also noticeable that there is no change in the channels as they pass across the transition from incised canyons (which Hobley et al.  showed were best described with a detachment-limited erosion model) to alluvial-fan sedimentation. In contrast to the batholith catchments, there is a high degree of variability in the molasse catchments, which also have higher Mχ values for opposing catchments. There are also changes associated with mapped structures within the Indus molasse such as the Choksti thrust (Sinclair and Jaffey, 2001). It is unlikely that the variation in Mχ values in the molasse has been caused by variations in base level along the Indus, because if this were the case, the variability in Mχ would be mirrored in the Ladakh batholith. Variability is more likely caused by changes in channel erodibility or changes in local erosion rates. Changing drainage areas due to divide migration can also alter Mχ values, but changing drainage area is unlikely to cause discontinuities in middle reaches of channels such as those seen in Figure 9D. Erodibility could be driving the difference in channel steepness, but 10Be erosion rates are significantly higher in the molasse, suggesting that differential erosion plays a role in setting the differences in channel steepness across the two geologic units.
We then turn our attention to possible structures (i.e., the Stock thrust) within the molasse that might lead to either increased local relative uplift (i.e., increased uplift relative to the Ladakh batholith) or changes in drainage area. If the molasse is being thrust toward the northeast, leading to motion of the drainage divide relative to the Indus, it would truncate the drainage area at the base of the catchment at the point of the thrust fault but would not affect the drainage area upstream. This is because the entire catchment would be advected to the north. On the other hand, if there were internal deformation within the molasse, in which drainage areas were systematically declining within the molasse, then according to Equation 1, χ would increase while elevation remained relatively constant, leading to a decrease in Mχ, which is the opposite of what we observe. A vertical component of thrusting would lead to increases in channel gradients and erosion rates across any faults, which is consistent with our observation of greater Mχ values in the molasse. This is corroborated by data from cosmogenic 10Be (see later section on cosmogenic nuclides).
We therefore find the most likely interpretation of the contrasting Mχ values between the molasse and the Ladakh batholith is the presence of at least one active thrust fault within the molasse. We cannot rule out the possibility that erodibility differences have caused the observed patterns, but given that erosion rates and fission-track data appear to correlate with the patterns in channel steepness (Fig. 9), we believe thrusting is the most likely explanation. We propose that the northeastward-vergent Stok thrust, as identified in the deformed terraces (Fig. 6), can be traced as an active structure along the range front at the head of the large alluvial fans that feed the Indus Valley (Fig. 8B), and that there is likely to be additional active displacement across other structures within the Indus molasse such as the Choksti thrust (van-Haver, 1984; Sinclair and Jaffey, 2001).
Sedimentary Evidence of Northward Migration of the Indus River Channel
In addition to the deformed terraces and steepened river profiles, there is sedimentary evidence to indicate that the course of the main Indus River channel has migrated northeastward through time.
The dissection of the T1 fan surface, described previously (Fig. 5B), exposes the internal stratigraphy of the fan, which reveals a unit composed of coarse boulder conglomerates with very well-rounded clasts up to 1.5 m in diameter, including multiple lithologies but with granodiorite from the Ladakh Range being dominant. Boulders and pebbles show strong imbrication indicating flow toward the northwest (i.e., parallel with and downstream from the modern Indus River). Exposures of this boulder conglomerate are seen in isolated locations higher up the fan, ∼1.2 km from the modern Indus River channel and 120 m higher (Fig. 5C).
Overlying the boulder conglomerate, there is a poorly structured gravel that includes angular clasts of the Indus molasse. These gravels are very poorly sorted, with some clasts greater than 1 m. The vague bedding dips gently down the direction of the dissected fan surface. In the middle of these gravels, there is a light cream-colored bedded and laminated siltstone, with interbeds of the gravels.
This upper succession represents deposits of the ancient alluvial fan interbedded with lake sediments that are traceable into the deformed T1 lake sediments described previously ∼3.7 km west-northwest from this location as unit 4 (Figs. 6 and 7). Underlying the alluvial gravel, the boulder conglomerates must represent the course of the Indus paleochannel prior to ca. 51 ka (oldest age of the T1 terrace from Phartiyal et al., 2005). The implication is that the modern Indus River channel has migrated northeastward since ca. 51 ka.
EROSION RATES ACROSS INDUS VALLEY
Having calculated rates of structural displacement across the Stok thrust, we synthesized published erosion rates from the upthrown side of the fault in order to evaluate the balance between vertical displacement rates and erosion rates. In addition, we compared these rates to the time-equivalent erosion on the opposite side of the Indus Valley from the Ladakh batholith, since the river morphologies suggest lower erosion rates. Published data on bedrock thermochronology and detrital cosmogenic 10Be are presented as a record of long-term (>106 yr) and short-term (<105 yr) data.
Thermochronology data on the cooling histories of rock samples within the top few kilometers of Earth’s surface in most mountain ranges can be used as an approximation of erosion rates (Reiners and Brandon, 2006). Apatite fission-track and apatite and zircon U-Th/He data have been extensively published from across the Ladakh region, with the majority of the analyses on the Ladakh batholith (e.g., Kirstein et al., 2006, 2009). We integrated these data with published values from the Indus molasse in the Stok Range (Clift et al., 2002; Sharma and Choubey, 1983).
The age-elevation data from the center of the Ladakh batholith as recorded by all three thermochronometers indicate rapid cooling at around 20 Ma (Kirstein et al., 2006), possibly linked to southward-vergent thrusting of the batholith. However, the lower-elevation interfluve promontories on the southwestern margin of the batholith nearest to the modern Indus River have older ages (Fig. 9). For example, the apatite fission-track ages range between ca. 35 and 30 Ma (Kirstein et al., 2009); this increase in age at lower elevations on the southern margin of the batholith remains to be explained.
Published apatite fission-track ages from the Indus molasse in the Zanskar Gorge record central ages of 13.7 ± 3.2 Ma and 13.8 ± 1.9 Ma (Clift et al., 2002). An additional age from further east was reported as being between 7 and 9 Ma (Sharma and Choubey, 1983). Assuming similar geothermal gradients across the Indus Valley, the contrast in ages (Fig. 9B) from the Indus molasse (ca. 14 Ma) to the southwestern margin of the Ladakh batholith (ca. 30–35 Ma) and into the core of the batholith (ca. 20 Ma) implies that the long-term erosion rates in the Indus molasse have been at least twice as fast relative to those in the batholith since at least ca. 14 Ma. The absolute erosion rates are difficult to assess due to lack of multiple thermochronometers and vertical profiles, but assuming a geothermal gradient of ∼30°/km, and closure temperature of 110 °C (e.g., Reiners and Brandon, 2006), then the likely erosion rates have been of the order of 0.1–0.3 m/k.y. We interpret the higher longer-term erosion rates in the Indus molasse to have resulted from Miocene to recent deformation, and the development of the Stok Range (Sinclair and Jaffey, 2001).
Concentrations of cosmogenically induced radionuclides such as 10Be and 26Al in quartz are routinely used for dating the period of exposure of a rock at the surface (Lal, 1991). Applications include dating boulders on glacial moraines (e.g., Brown et al., 1991) and fluvial bedrock strath terraces (e.g., Burbank et al., 1996b). Additionally, cosmogenic radionuclides measured from quartz sand in river catchments can be used to estimate averaged catchment-wide erosion rates (Lal and Arnold, 1985; Brown et al., 1995). This method fills the gap between traditional erosion estimates determined from measured river sediment loads (Schaller et al., 2001) and long time scales approximated using thermochronology.
Analysis of cosmogenic 10Be from sediment across the Ladakh batholith has demonstrated catchment-averaged erosion rates ranging from ∼0.02 to 0.08 m k.y.–1, with ∼10% error (Dortch et al., 2011b). These rates are slow compared to the mean for the Himalaya mountain range as a whole, which is ∼1.0 m k.y.–1 (Lupker et al., 2012). The catchments along the southern side of the Ladakh batholith have calculated mean averaged erosion rates of between ∼0.02 m k.y.–1 and 0.04 m k.y.–1 (Fig. 10; Dortch et al., 2011b; Munack et al., 2014). Further detrital 10Be data from the Stok Range record rates of 0.07–0.09 m k.y.–1 (Munack et al., 2014); this supports the fission-track thermochronology by indicating erosion rates of the Indus molasse that are approximately twice as fast as those over the batholith on the opposing side of the Indus Valley (Fig. 9B).
INTERPRETATION OF RESULTS (FIG.11)
These results confirm that structural shortening is taking place across the present Indus Valley, with horizontal displacement rates of at least 0.21 m k.y.–1, which represent just a small fraction (<2%) of the total shortening across the Himalaya in this region. An assumption in this calculation is that the recurrence interval between episodes of fault slip on the Stok thrust are approximately equal to the period since the last displacement. If the recurrence intervals are longer, then the time-averaged slip rate on the thrust would be slower.
The deformation of the Stok Range has resulted in a steepening of river channels with higher erosion rates and a higher sediment flux into the Indus Valley relative to the opposing tributaries that drain the Ladakh batholith. The steepened river channels enable the erosional response to variable rock uplift relative to the Indus River valley to be mapped to the east of the observed thrust, indicating that the Stok thrust has been active along the northeastern margin of the mountain front.
Whether the horizontal displacement rate across the Indus Valley is sufficient to permanently offset the drainage divides depends on the ability of the erosive processes to counter the topographic displacement induced by the deformation. In bedrock channel networks, erosion is enhanced by the propagation of knick zones (or knickpoints) up to the head of the catchments and the consequent response of hillslopes. In geometric terms, a river catchment’s ability to fully recover its form during shortening across a thrust fault can be simplified to a ratio of the vertical rock displacement rate (Vv; this being displacement relative to the footwall block) versus the erosion rate in the hanging wall of the thrust (Fig. 1). The vertical rock displacement rate is a function of the horizontal displacement rate (Vh) multiplied by the combined tangents of the dip of the thrust (β) and the mean topographic slope (α) from ridgeline to tributary outlet (Fig. 1). The Stok thrust has a measured dip at the surface of 37°, and the mean surface slope of the Stok Range from ridge crest to valley floor is ∼8° (Fig. 9A). Therefore, a horizontal displacement rate of 0.21 m k.y.–1 equates to a vertical displacement rate of ∼0.19 m k.y.–1.
The catchments that drain the hanging wall of the Stok thrust are sourced from the Indus molasse, where the erosion rates measured from 10Be concentrations range from ∼0.07 to 0.09 m k.y.–1. This implies that more than half of the vertical component of displacement, and proportionately half of the horizontal component of displacement, on the fault is converted into topographic displacement at the surface, with the rest being eroded. The implication is that the drainage divide that forms the spine of the Stok Range must be migrating toward the Indus River at approximately half the rate of horizontal displacement on the Stok thrust, equating to ∼0.1 m k.y.–1. However, given the presence of knick zones up the Stok Range catchments (Fig. 9C), the 10Be concentrations are likely recording a mixture of higher erosion rates (lower 10Be concentrations) below the knick zones and lower rates (higher 10Be concentrations) above the knick zones, where the kinematic wave of accelerated incision has not yet reached. For this additional reason, it is possible that the calculation for divide migration is an underestimate, and that the true value is likely to lie somewhere between 0.1 m k.y.–1 and the horizontal fault displacement rate of 0.21 m k.y.–1. If the knickpoints were able to propagate to the drainage divide prior to subsequent shortening on the fault, then the averaged erosion rates for the catchments would be considerably higher, and they could fully reequilibrate the channel profile and so counter the fault displacement rates. However, the presence of multiple channel segments of variable steepness values (Fig. 9D) bounded locally by knickpoints (Fig. 9C) suggests reequilibration of the channels is not achieved between episodes of shortening across the fault.
Greater fault displacement rates compared to erosion rates in the hanging wall of the Stok thrust demonstrate that this topographic form is evolving, and that the elevation contrast from outlet to drainage divide across the Stok Range (catchment relief) is likely to be increasing over long time scales. If we assume the elevation of the Indus Valley is constant, then it would suggest that catchment relief is growing at a similar rate to the divide migration rate, i.e., ∼0.1 m k.y.–1. However, the elevation of the Indus Valley relative to the surrounding tributaries has fluctuated, as recorded in the documented alluvial terraces, but the present elevation of the Indus River is up to 100 m lower than it was ca. 50 ka (Fig. 4B). Based on this evidence, it is hard to conclude whether the long-term elevation of the Indus River channel is rising or falling relative to the deforming Indus molasse of the Stok Range.
As sediment flux increases with relief and channel steepening, so must have the alluvial fans that drain into the Indus Valley off the Stok Range expanded. The expansion of alluvial fans from this side of the valley has forced the present-day Indus channel to migrate laterally toward the opposing valley margin against the rock promontories of the batholith (Fig. 11). This interpretation is supported by the presence of Indus River boulder conglomerates exposed beneath the present fans that drain the Indus molasse (Figs. 5B and 5C). Another consequence of the asymmetry in erosion and sediment flux across the Indus Valley is the aggradation of alluvial fans within the valleys of the Ladakh batholith. This aggradation has resulted in some isolated hills or inselbergs of granodiorite that once formed parts of interfluve ridges, but are now buried in alluvium and topographically detached from the range.
While this study has focused on the tectonic driver for topographic narrowing of the Indus Valley, it is clear that asymmetry of erosion rates driven by lithology and climate will also influence divide migration. In the case of the Indus Valley, the divide that runs along the Ladakh batholith has also experienced a strong asymmetry in glacial erosion, with headwall retreat rates of 0.18–0.6 m k.y.–1 in the northeastward-facing glaciated catchments (Dortch et al., 2011a). A similar asymmetry exists across the drainage divide of the Stok Range, with glaciers only present on the northeast-facing catchments; however, the rates of headwall retreat have not been measured in these catchments. A working assumption is that the rates of headwall retreat across these two drainage divides is comparable. Dating the onset of glaciations in these settings has not been possible, but 10Be exposure ages extend back to ca. 300 ka in the region (Dortch et al., 2013), and it is probable that glaciers have existed in this setting for a great deal longer. Therefore, the tectonic convergence of drainage divides must have been modified by the glacial signal. In the case of the Indus Valley, if the two divides on the batholith and Stok Range are equally influenced with the same sense of motion, then the impact on the relative convergence between the drainage divides is zero. However, it is recognized that this remains a component of significant uncertainty.
An additional consequence of the thrust deformation driving topographic shortening across the valley is the increased susceptibility to damming of the valley (e.g., Bürgisser et al., 1982; Blöthe et al., 2014). We have been able to demonstrate that the thickest lacustrine terrace in the Leh Valley was most likely caused by thrust motion at ca. 38 ka on the Stok thrust, leading to debris flows and consequent damming of the valley. The younger T2 terrace can also be correlated downstream to similar age deposits that incorporate numerous mass-flow deposits (e.g., Blöthe et al., 2014).
Here, we have documented the structural, topographic, and surface process response to slow horizontal displacements across a single valley. We would expect similar processes to take place simultaneously across numerous structurally defined, strike-parallel (longitudinal) valleys in any large mountain range. In the Himalaya, rivers such as the Tsangpo and Shyok, and the upper reaches of the Kosi, Sutlej, and Karnali all run parallel to structures that may be actively modifying catchment form in a similar way to the Indus case presented here. This is particularly relevant where active shortening occurs in regions of relatively low erosion rates, such as in the lee of the Himalaya.
(1) Through OSL dating and analysis of Quaternary terraces in the Indus Valley, Ladakh, we demonstrated that the southwestwardly dipping Stok thrust, which represents a lateral continuation of the Great Counter thrust in the Himalaya, has been active from ca. 42 ka, resulting in ∼10 m of shortening. The displacement on this structure resulted in debris flows blocking the valley, and the formation of a lake in this part of the Indus Valley.
(2) Mapping of river channel steepness using the chi parameter in the hanging wall of the Stok thrust indicates that its recent activity is laterally traceable at least 80 km southeastward along the valley. The variably steepened channels of the Stok Range contrast with the relatively steady, lower-gradient (normalized for area) channels that drain the Ladakh batholith on the opposing side of the Indus Valley.
(3) Erosion rates as recorded by low-temperature thermochronology and detrital 10Be concentrations from river sediment are approximately twice as fast over the Stok Range (up to 0.09 m k.y.–1) relative to the Ladakh batholith (up to 0.04 m k.y.–1).
(4) Deposits of the Indus River channel buried beneath the alluvial fans sourced from the Indus molasse testify to the northeastward migration of the present river channel. The interpreted mechanism is that relatively high sediment yield from the Stok Range has forced the course of the present channel against the opposing valley margin formed by the batholith. In addition, the high sediment yield has forced fluvial base levels to rise over the batholith, blanketing the interfluve bedrock ridges with alluvium.
(5) The contrast in the vertical rock displacement rate of the hanging wall of the Stok thrust (∼0.19 m k.y.–1) versus the erosion rate (<0.09 m k.y.–1) requires a change in surface topography. The fact that approximately half of the vertical rock displacement is countered by erosion implies that approximately half of the structural shortening is recorded as topographic convergence of drainage divides across the valley (assuming glacial headwall retreat is comparable on both drainage divides). This recorded deformation of the Indus River valley at rates of ∼0.1 m k.y.–1 represents the first documentation of the processes and consequent topographic and sedimentological record of narrowing of a major longitudinal river valley in an active mountain range.
Funding for this project was provided from the Natural Environment Research Council (NERC grants NE/I017747/1 and NE/J012750/1), and the Carnegie Trust for the Universities of Scotland. The authors are thankful to Fida Mitoo Hussein for field logistics. Fida has sadly since died, and we would like to acknowledge his help for much of the work in this region. We thank Mark Brandon, Frank Pazzaglia and an anonymous reviewer for valuable comments on the initial submission.
- Received 2 October 2015.
- Revision received 11 March 2016.
- Accepted 26 July 2016.
- © 2016 Geological Society of America