Segmented strain accumulation in the High Himalaya expressed in river channel steepness Open Access

The Himalayan thrust wedge is 2500 km long, 250 km wide, and up to 50 km thick below southern Tibet (Nábělek et al., 2009) (Fig. 1). The wedge is floored by an intracontinental thrust fault separating India from Asia known as the Main Himalayan thrust (MHT), whose behavior gradually changes from an aseismically sliding décollement beneath southern Tibet to the Himalayan seismogenic zone across a 50–100-km-wide brittle-plastic transition (BPT) beneath the High Himalaya (Bilham et al., 2001; Ader et al., 2012; Cannon and Murphy, 2014; Stevens and Avouac, 2015) (Fig. 2). The BPT corresponds to the boundary between the strongly coupled portion of the MHT, which only moves during earthquakes, and the weakly coupled portion (Stevens and Avouac, 2015). Clusters of steep river channels are located at the BPT (Fig. 2). Clusters are offset from each other across strain partitioning structures that have been independently proposed to separate the range into orthogonally convergent and obliquely convergent sectors (Drukpa et al., 2006; Velasco et al., 2007; Styron et al., 2011; Murphy et al., 2014) (Fig. 3).

A record of Himalayan strain accumulation is encoded in the landscape through a competition between erosion and rock uplift, and can be visualized by mapping normalized river channel steepness (k_sn), a proxy for rock uplift over 10⁵ yr. Currently, high topography, high-k_sn rivers, rapid denudation rates (Godard et al., 2014; Morell et al., 2017), rapid interseismic strain rates (from GPS velocities) (Ader et al., 2012; Bilham et al., 2001; Stevens and Avouac, 2015), orogen-parallel belts of seismicity (Bollinger et al., 2007; Pandey et al., 1995, 1999), and MHT thrust-ramp duplexes are all coincident with the inferred MHT BPT (Fig. 1B) (Pandey et al., 1999; Bollinger et al., 2004a; Herman et al. 2010; Cannon and Murphy, 2014; McQuarrie et al., 2014; Morell et al., 2015; Harvey et al., 2015) tracing the northern limit of the strongly coupled plate boundary (Stevens and Avouac, 2015) (Fig. 2). This suggests that the position of high topography within the thrust wedge is invariant over 10⁵ yr time scales, and that it is a function of thrust wedge geometry, plate boundary thermal structure, and stored elastic strain. As described below, the strongest correlations found in our analysis are between k_sn, MHT coupling (–0.6), and earthquake density (0.6) (excluding regions with too few modern earthquakes for statistical significance), and topography (0.6) highlighting that in active orogens, k_sn records strain accumulation (modified by substrate erodibility, channel width, sediment load, and climate). All other parameters are correlated with k_sn.

The MHT’s surface trace is manifest as a series of thrust-related folds and faults at the toe of the wedge, which form the first spatial, southernmost topographic expression of the Himalaya and are collectively known as the Main Frontal thrust (Fig. 1). The rate of thrusting across the Main Frontal thrust in central Nepal, documented by offset and folded river terraces (Lavé and Avouac, 2001), is within uncertainty of the rate of aseismic sliding recorded by GPS stations in south-central Tibet (Ader et al., 2012; Bilham et al., 2001, 1997; Jackson and Bilham, 1994). This has been interpreted by some researchers to signify that the wedge does not deform internally, implying that rock uplift within the wedge is a passive product of fluxing material over MHT ramps (Cattin and Avouac, 2000; Pandey et al., 1995). However, the suggestion of an entirely passive wedge is at odds with other observations. Recent work has highlighted that the A.D. 1505 earthquake in west Nepal may have ruptured multiple splays within the wedge (West Nepal fault system, Main Boundary thrust) (Hossler et al., 2016; Murphy et al., 2014; Harvey et al., 2015), while aftershocks of the A.D. 2015 Gorkha earthquake appear to have ruptured internal thrust wedge faults at the base of the High Himalaya (Whipple et al., 2016; Hoste-Colomer et al., 2016). Quaternary thrust faulting at the base of the High Himalaya has been proposed on the basis of field relationships and thermochronology (Hodges et al., 2004; Wobus et al., 2005), although evidence of recent surface rupturing thrusts at the base of the High Himalaya has been elusive. Our k_sn map, as well as thermokinematic models of late Miocene to recent times and structural restorations of the frontal portion of the wedge are consistent with underplating of Indian crust along a foreland-migrating ramp (Bollinger et al., 2004a; Herman et al., 2010; McQuarrie and Ehlers, 2015; Robinson, 2008).

GPS measurements show that the most-rapid interseismic strain accumulation is located near or above (north of) the 3500 m topographic contour (Ader et al., 2012; Bilham et al., 2001), which has been used to model MHT coupling (Stevens and Avouac, 2015). We use this model to examine the relationship between river channel steepness and interseismic strain (Fig. 2). There is a nonlinear relationship between k_sn and rock-uplift rate (Kirby and Whipple, 2001, 2012), where the time required to impart a k_sn signal to the landscape is dependent on the long-term (multiple seismic cycles) rock-uplift rate. In the Himalaya, vertical interseismic strain rates are as much as 4–7 mm yr^–1 (Jackson and Bilham, 1994; Ader et al., 2012; Grandin et al., 2012; Stevens and Avouac, 2015), however the 25 April 2015 Gorkha earthquake demonstrated that much of this interseismic uplift is elastic and destined to be transferred into horizontal and vertical motion of the strongly coupled portion of the range during plate boundary earthquakes (Elliott et al., 2016; Lindsey et al., 2015).

During the Gorkha earthquake there was a ∼0.5 m lowering of the High Himalaya and 1 m of uplift for the Lesser Himalaya (Elliott et al., 2016; Lindsey et al., 2015), showing that the long-term uplift rate of the High Himalaya may be a fraction of that measured by GPS, interferometric synthetic-aperture radar (InSAR), and spirit leveling. The map presented here documents clusters of high-k_sn rivers, most of which straddle the boundary between the strongly coupled and weakly coupled portions of the MHT (Stevens and Avouac, 2015). In regions with a high rate of seismicity, clusters of earthquake epicenters are located south of (foreland to) the steepest river channels (Fig. 2). We follow Cattin and Avouac (2000) and Bollinger et al. (2004b) in interpreting peak seismicity to reflect the leading edge of the BPT and, based on our data, suggest that peak k_sn reflects focused elastic strain in the upper crust and plastic deformation of the lower to middle crust (Cannon and Murphy, 2014). Here within the BPT, slip is fed from the plastically deforming lower crust beneath southern Tibet to MHT thrust-ramp duplexes in the mid- to upper crust via underplating or tectonic wedging. The MHT BPT is a mechanical boundary whose location within the thrust wedge is a function temperature along the subduction interface (MHT). The southernmost limit of the BPT, closest to the toe of the thrust wedge, begins where the MHT encounters the 350 °C isotherm (Herman et al., 2010), at which point quartz begins to deform plastically (Hirth and Tullis, 1992). The percentage of deformation accommodated by crystal-plastic creep increases with temperature as the dominant mineral constituents begin to deform plastically (feldspar begins to deform ductily at 600 °C) (Tullis and Yund, 1992). The width of the MHT BPT depends on the geothermal gradient and the local dip of the MHT. The thrust wedge north of the BPT acts as a dashpot (a device which resists motion via viscous friction), storing strain until it overcomes static friction on the strongly coupled portion of the MHT, generating a great earthquake (Bilham et al., 2001; Avouac, 2003; Feldl and Bilham, 2006; Cannon and Murphy, 2014). Either some amount of this strain is permanent or the high-k_sn signal present across the High Himalaya represents a purely elastic bulge maintained by ongoing India-Asia convergence as suggested by Meade (2010). If true, then the highest topography of the Himalaya is dynamically maintained by plate convergence. While this idea is unconventional, an elastic strain reservoir equivalent to four M_w 8 earthquakes or a single M_w 9 earthquake has been inferred to exist across the range by comparison of the seismic energy released from all known great Himalayan earthquakes of the past 500 yr with convergence over the same time interval (Bilham and Ambraseys, 2005; Stevens and Avouac, 2016). While the percentage of accumulated strain released in plate boundary earthquakes versus absorbed by internal shortening structures remains uncertain (Bilham and Ambraseys, 2005; Cannon and Murphy, 2014; Meade, 2010; Stevens and Avouac, 2015; Taylor, 2016), the identical rate of frontal anticline growth and convergence of southern Tibet and India (Lavé and Avouac, 2001) indicates that the amount of strain available for internal deformation is small.

Bilham et al. (2001) estimated that up to 10% of the strain accumulated in the High Himalaya is available for crustal thickening. This slow uplift of the High Himalaya is evident in the k_sn map, as clusters of high-k_sn rivers in the High Himalaya are offset from each other across Himalayan-Tibetan rifts. Himalayan river channel steepness can be used as a proxy for strain accumulation and provides an intermediary time step between the decadal GPS record, the roughly 1000 yr historical record of great earthquakes, and the ±0.5 m.y. record of low-temperature thermochronology.

The surface bedrock geology of the High Himalaya is largely the result of emplacement of the metamorphic core of the orogen between the Main Central thrust and the South Tibetan detachment (Fig. 1) during the early Miocene (Searle and Godin, 2003; Edwards et al., 1999; Murphy and Harrison, 1999; Vannay et al., 2004; Thiede et al., 2006; Kellett et al., 2010). In central Nepal, these structures appear to maintain a low degree of activity today (Hodges et al., 2004; McDermott et al., 2013; Wobus et al., 2005), however no recent surface ruptures have been found, and cooling histories and river channel steepness data sets have been variously interpreted to represent out-of-sequence faulting, duplexing, or passive uplift over MHT ramps (Pandey et al., 1995; Cattin and Avouac, 2000; Avouac, 2003; Bollinger et al., 2004a; Wobus et al., 2005; Herman et al., 2010; McDermott et al., 2013, 2015; Hubbard et al., 2016).

In west Nepal, the South Tibetan detachment is folded (Fig. 1) and defines the Dolpo anticline, a regional structure 400 km long with at least 9 km of structural relief whose frontal limb coincides with a cluster of high-k_sn channels interpreted to be within the BPT in west Nepal (Fuchs, 1964; Cannon and Murphy, 2014). Folding initiated after the cessation of South Tibetan detachment motion, regionally constrained to 16–19 Ma (Searle and Godin, 2003; Godin et al., 2006b). The concentration of high-k_sn rivers and young zircon (U-Th)/He ages (8.0 ± 1.3 to 2.6 ± 0.7 Ma) (McCallister et al., 2014) at Gurla Mandhata (a core complex marking the western limit of the Dolpo anticline) (Murphy et al., 2002; Murphy and Copeland, 2005) suggests that hinterland folding initiated in the late Miocene and is ongoing today.

In Bhutan, the South Tibetan detachment was active 23–16 Ma (Kellett et al., 2010) or 20.5–14 Ma (Cooper et al., 2015) and appears to be folded into a broad anticline that spans the northern provinces of Punakha, Bumthang, and Kurtoed (Grujic et al., 2002; Long et al., 2011; McQuarrie et al., 2014). The anticline is collocated with the BPT interpreted from our k_sn data. Low-temperature thermochronology of the northern Bhutan Himalaya is largely lacking from the literature, but Adams et al. (2013) presented apatite (U-Th)/He ages (2.4 ± 0.4 to 3.6 ± 1.2 Ma) from the southern edge of the northern Bhutan high-k_sn cluster, which reflects recent exhumation of the Bhutanese hinterland. The anticlines overlying the BPT of west Nepal and Bhutan, interpreted from k_sn data, record recent exhumation, and are bound on one side by major extensional structures, the Gurla Mandhata core complex and the Chomolari detachment (Yadong-Gulu rift) (Murphy et al., 2002; Wu et al., 1998). In both cases, folding appears synchronous with the onset of extension, suggesting that the initiation of hinterland thrust duplexes is kinematically linked with thickening of the lower crust in a constrictional strain field (Murphy and Copeland, 2005; Sundell et al., 2013).

Normalized river-channel steepness (k_sn) analysis was conducted using the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model V002 with 30 m resolution from 78°E to 92°E. Although the ASTER digital elevation model has better resolution than previous models, it contains numerous holes and data issues. Deeply incised river canyons sometimes generate “digital dams”, which artificially elevate the river channel. Reaches of rivers with bad data are recognizable on longitudinal profiles (plots of elevation versus downstream distance) by their jagged stair-step appearance. Good k_sn fits are not possible for these reaches, hence they were not analyzed. For this reason, deep canyons on along-trunk channels within the High Himalaya were commonly not able to be analyzed. To complete our k_sn map in these regions, we analyzed tributary channels instead of trunk channels. This is the advantage of using normalized channel steepness, as it allows the user to quantitatively compare channels with different drainage areas.

(Kirby and Whipple, 2001; Wobus et al., 2006), in which S is slope, A is upstream drainage area, k_sn is normalized channel steepness, and θ_ref is reference concavity (0.45). The digital elevation model was hydrologically corrected in ArcGIS software (https://www.arcgis.com/) using the “fill” function in the hydrology toolbox to fill “digital pits” in the data. Then a flow-direction grid was created by determining the downslope direction at each cell using the “flow direction” function. Finally, using the flow-direction grid as an input, an accumulation grid was created by calculating the accumulated flow at each cell in the digital elevation model using the “flow accumulation” function. Once these hydrologic corrections were completed, the analysis could proceed. In ArcGIS, analysis parameters were selected for each session. In areas where the data quality is good, a 250 m smoothing window was used; in areas with poor data quality, a 500 m smoothing window was used. A 0.45 reference concavity was used for all sessions. Site-specific reference concavity can be calculated, but studies examining its effect on the analysis have all found that the choice of reference concavity changes the values obtained, but not the pattern of high and low channel steepness (Kirby et al., 2003; Whipple, 2004; Tyagi et al., 2009; Olen et al., 2015). Studies of k_sn in mountainous regions use reference concavities between 0.4 and 0.5, with 0.45 being the most common choice (Kirby and Whipple, 2001, 2012; Kirby et al., 2003; Wobus et al., 2003, 2006; Duvall et al., 2004; Hodges et al., 2004; Korup, 2006; Clark and Bilham, 2008; Robl et al., 2008; Cyr et al., 2010; Allen et al., 2013; Adams et al., 2016; Harvey et al., 2015). After the parameters were set, a point on a river near its headwaters was selected, and its channel data (elevation, drainage area, and gradient at each point along the channel) were exported to Matlab software (https://www.mathworks.com/products/matlab.html). In Matlab, reaches of rivers with uniform gradient were visually selected on the longitudinal profile, exported to ArcGIS, and plotted on the digital elevation model. Deviations from steady-state longitudinal profiles can be produced by changes in rock-uplift rate, changes in substrate erodibility, and climatic gradients, so the results of the analysis must be evaluated for each of these forcing mechanisms.

Strong range-normal climatic gradients exist across the Himalaya, where the Lesser Himalaya receives two to four times more annual rainfall than the High Himalaya. Significant along-strike climatic gradients also exist across the arc, where the eastern Himalaya is significantly wetter than the western Himalaya. Despite these differences, the High Himalaya can be generalized as being contained within a single relatively dry climate and the Lesser Himalaya as belonging to a wetter climate. Because the emphasis of this study is strain accumulation in the High Himalaya, we have not weighted the flow-accumulation grid for climatic variations. Doing so would lower k_sn values slightly in drier regions and raise them slightly in wetter regions, but would not change the overall pattern. Determining appropriate weighting factors for the various arc-normal and arc-parallel climatic variations present across the Himalaya is outside the scope of the current study, and given the low level of correlation between precipitation rates and river channel steepness, it would be unlikely to significantly affect the results of the analysis.

The relationship between rock uplift, upstream drainage area, and river channel morphology in glaciated regions cannot be effectively modeled using the k_sn index. Measurements of modern Himalayan glacial equilibrium-line altitudes are broadly bracketed at 4000–5000 m (Owen and Benn, 2005). The Himalaya has a long complicated history of glaciations during which equilibrium-line altitudes have fluctuated considerably. During the last glacial maximum, many Himalayan glaciers reached below 1000 m elevation, and the toe of some nearly reached the Gangetic plain (Owen and Benn, 2005, and references therein). For this reason we, like in all other studies of Himalayan river channel steepness, cannot exclude all landscapes that have had any glacial influence. Instead we tried to exclude active glacial landscapes, while not eliminating more of the High Himalaya than necessary, by clipping all raster data sets excluding elevations above 4500 m. An unavoidable consequence of including areas that experienced past glaciation is the apparent oversteepening of the lowermost portion of tributaries where they join a formerly glaciated trunk channel. Widening of the trunk channel during glaciation produced a knickpoint in the tributary channel at their junction. This is seen in our data as a short (<1 km) anomalously steep segment at many high-elevation trunk-tributary junctions. However, the effect of these glacial artifacts on the orogen-scale map is small due to their limited spatial extent.

To test the dependency of k_sn on interseismic versus coseismic uplift and how it relates to the generation of relief, we calculated correlation coefficients for k_sn versus MHT coupling, earthquake epicenter density, and topography. We evaluated the contributions of lithology (proxy for substrate erodibility) and precipitation by generating correlation coefficients for k_sn versus our geologic map [compiled from DeCelles et al. (1998); Long et al. (2011); Mitra et al. (2010); Yin (2006); Webb et al. (2007); and He et al. (2015)] and a 50 yr interpolation of precipitation data (Hijmans et al., 2005) (Supplemental Figs. S1–S3¹). The statistical analysis was conducted using ArcGIS and Microsoft Excel software. The k_sn data constitute a vector data set; comparing it quantitatively with raster data sets like MHT coupling, topography, and precipitation required converting the individual vector shapefiles to a single raster file. We used the natural neighbors interpolation algorithm with a 1 km grid for the conversion. The final step in preparing the k_sn data for display was to contour the data to highlight continuous clusters of high-k_sn rivers. We define strain accumulation segments as regions near the BPT whose k_sn values are consistently >200 (Fig. 2).

Similarly, earthquake epicenters of M_w >2.5 from A.D. 1966 to 2016 from the U.S. Geological Survey National Earthquake Information Center (NEIC) database (https://earthquake.usgs.gov/earthquakes/search/) were converted to a raster grid using the point density algorithm with an output cell size of 0.05° and circular search radius of 0.5° (Supplemental Figs. S1–S3 [^{footnote 1}]). Data were then extracted from each data set along the same line using the “interpolate line” function on the “3D Analyst” toolbar in ArcGIS, followed by the “profile graph” function.

Hoek and Brown created their failure criterion as a way to estimate the strength of rock masses based on lithology through testing of uniaxial and tensile strength of different lithologies. We took the m_i values of all lithologies within each Himalayan lithotectonic unit and averaged them to reflect average rock mass strength for each (Table 1), and then normalized the values to 1 for comparison with the other data sets. The resulting erodibility index values are: 0.5 for the Siwaliks, 0.7 for the Lesser Himalaya sequence, 1.0 for the Greater Himalaya sequence, and 0.5 for the Tethyan sedimentary sequence, which essentially separates the Himalaya into soft, intermediate, and hard rocks. To characterize the relationship between k_sn, MHT coupling, substrate erodibility, and precipitation, we analyzed four orogen-normal profiles spaced at 25 km intervals for ten regions spanning the central Himalaya from Uttarakhand to Bhutan (Supplemental Figs. S1–S3 [^{footnote 1}]). The results reported for the entire range (Table 2) are an average of all ten regions.

The results of the k_sn analysis and the statistical analysis of potential forcing mechanisms and their effects are presented in the following sections. First, the role of precipitation and lithology in determining river channel steepness is evaluated. Next, the contribution of strain accumulation (interseismic and coseismic), in the form of MHT coupling and earthquake epicenter density, is examined. Finally, the relationship between topography, k_sn, earthquakes, and drainage networks is investigated.

A map of average rainfall interpolated for the past 50 yr (Hijmans et al., 2005) was used to generate correlation coefficients for k_sn and precipitation. The average r (Pearson correlation coefficient) for the entire Himalaya is –0.3 with a p < 0.005 (p-values measure the likelihood that the correlation results from the null hypothesis; any p < 0.05 is considered statistically robust) (Fig. 4; Supplemental Figs. S1–S3 [^{footnote 1}]). No relationship was found between the average rate of rainfall and the strength of the correlation. The low degree of average correlation reveals that either the spatial pattern of Himalayan precipitation has undergone significant changes over the last ten thousand to hundred thousand years, or that erosion rates are decoupled from precipitation rates as suggested by recent studies of the Alpine-Himalayan orogen (Burbank et al., 2003; Godard et al., 2014; Adams et al., 2015, 2016; Forte et al., 2016).

This analysis reveals a moderate degree of positive correlation (average r = 0.5, p < 0.05) for k_sn and lithology (Fig. 4; Supplemental Figs. S1–S3 [^{footnote 1}]). Regionally, the correlation is strong (0.6–0.8) across east Nepal and Bhutan (see Table 2), demonstrating the effect that the juxtaposition of lithologies with different erodibilities has on river channel steepness. The moderate degree of average correlation between k_sn and litholotectonic unit is evidence of the role substrate erodibility plays in modulating the k_sn signal; while it can be locally significant, e.g., capstone plateaus protecting more erodible substrates, regionally it appears secondary to rock-uplift rate in determining the pattern of river channel steepness.

The MHT coupling model of Stevens and Avouac (2015) is based primarily on GPS data but also uses the northern limit of clusters of earthquake epicenters to aid in locating where the plate boundary becomes uncoupled. MHT coupling displays a moderately strong correlation with k_sn (r = –0.6, p < 0.05) (Fig. 4; Supplemental Figs. S1–S3 [^{footnote 1}]). The correlation between the degree of MHT coupling and river channel steepness suggests that in the Himalaya, rock-uplift rate is the primary driver of river channel steepening, and relates it directly to interseismic strain. Tectonic steepening of rivers takes place over hundreds of thousands of years and is the net result of the interaction between interseismic, coseismic, and post-seismic strain accumulation and erosion; the relation between k_sn and coseismic strain release will be examined in the next section. Regionally, this correlation is strong across the western Himalaya (average r = –0.7) and becomes weak in the east (average r = –0.3) (see Table 2). Whether this is related to the paucity of GPS stations in easternmost Nepal, Sikkim (northeastern India), Bhutan, and Arunachal Pradesh (northeasternmost India) is unknown. While the GPS data upon which the MHT coupling model is based were recorded over tens of years, it takes far longer for the k_sn signal to be recorded in the landscape. While the exact amount of time required to impart a k_sn signal to the landscape is unknown, it is certainly over numerous seismic cycles. The strength of correlation between decadal interseismic elastic strain and k_sn demonstrates that strain accumulation has been localized in its current location for a period of time on the order of 100,000 yr. Our analysis suggests that of the three forcing mechanisms thought to influence river channel steepness (substrate erodibility, erosion, and rock-uplift rate), in the Himalaya, rock-uplift rate is the primary driver, with secondary influence from lithology, and to a more minor extent by the pattern of modern precipitation. Our correlations agree well with recent studies highlighting the dominance of tectonic processes over precipitation gradients in driving rock uplift (Burbank et al., 2003; Godard et al., 2014; Adams et al., 2015, 2016; Forte et al., 2016).

Earthquakes occur when accumulated interseismic strain overcomes static friction on the strongly coupled portion of a fault. Stored interseismic (elastic) strain is released, resulting in coseismic displacement of material and inelastic deformation taking the form of crustal thickening, thinning, or lateral offset. The spatial density of earthquake epicenters is used here as a representation of the amount and location of coseismic energy released. We use epicenters from the NEIC catalog from 1 January 1966 to 1 January 2016 because it is the longest continuous record available for the entire Himalayan arc. The 50 yr record of earthquakes used in our analysis is impacted by the presence of the 2015 Gorkha earthquake and its sequence of aftershocks, which dominate the east Nepal segment. Similarly, moderate-magnitude earthquakes and their aftershock sequences create local epicenter “hotspots” across the arc. Along-strike variations in the distribution of modern earthquakes should be viewed as a snapshot of coseismic energy release; we expect that given a sufficiently long record, earthquake clusters near the BPT will span the entire arc. We acknowledge that the 50 yr record used in this analysis is only a meaningful metric of coseismic energy release in regions where there have been modern earthquake sequences, and do not mean to imply that the lack of modern earthquakes in intervening regions suggests a fundamental difference in long-term strain accumulation. Instead we rely on the 10⁵ yr record of strain accumulation provided by k_sn mapping to denote regional differences. River channel steepness is impacted by both inter- and coseismic strain and as such should correlate well with any record of strain accumulation and release. This led us to test the relationship of k_sn with earthquake epicenter density. In regions with a sufficient number of earthquakes to produce a statistically significant result, epicenter density has a moderately strong correlation with k_sn (r = 0.6, p < 0.05) and relates river channel steepness directly to coseismic strain release and concomitant inelastic crustal deformation (Fig. 4; and Supplemental Figs. S1–S3 [^{footnote 1}]). Eight out of the ten regions analyzed met this criterion; the remaining two (Arun river area) and west Bhutan) were not used in calculating the average correlation. The moderately high degree of correlation between k_sn and earthquake epicenter density is demonstrative of the link between the steepening of river channels and crustal thickening by internal deformation of the high Himalaya. In regions with sufficient earthquakes for a robust statistical analysis, epicenter density peaks spatially southward of (foreland to) k_sn, supporting the idea that the High Himalaya acts as an elastic strain reservoir driving earthquakes on the strongly coupled portion of the plate boundary, and also demonstrates that rock uplift during earthquakes contributes to the k_sn signal.

Topography and k_sn display a moderately strong correlation (r = 0.6, p < 0.005), with the strongest correlation in the central and east Nepal segments (average r = 0.7) (Fig. 4; Supplemental Figs. S1–S3 [^{footnote 1}]). This is expected, as the same forces primarily responsible for steepening stream channels (stored elastic strain and crustal thickening) also generate the high topography. The segments with the strongest correlation between k_sn and topography are the same segments that host all of the Himalayan peaks >8000 m (excluding the western syntaxis), suggesting that these segments also accumulate strain more rapidly than the rest of the Himalayan range.

Over its 2500 km length, regions of rapid rock uplift in the High Himalaya are characterized by spatial clusters of high-k_sn rivers that are offset from one another across segment boundaries that coincide with linked strike-slip Himalayan-Tibetan rift systems. Nearly all high-k_sn clusters overlie weakly coupled portions of the MHT located 60–120 km north of the range front (Figs. 2 and 5). Only west Nepal and Bhutan host high-k_sn clusters that overlie strongly coupled portions of the MHT, 30–50 km from the range front (Figs. 2 and 5). The spatial dimensions of the high-k_sn clusters vary in width, length, and distance from the range front and are summarized in Table 3.

High-k_sn clusters overlying weakly coupled portions of the MHT are located at variable distances (60–120 km) from the range front, but nearly all are between 70 km and 100 km. Neighboring segments host clusters that are offset by 40–60 km. This indicates that whatever mechanism(s) (passive uplift above MHT ramps, out-of-sequence thrusting, duplexing, tectonic wedging) are responsible for rapid rock uplift within the clusters, it is not purely a function of distance from the range front as would be expected if the geometry of the MHT was uniform along strike. The segment boundaries are collocated with the Himalayan-Tibetan rift systems and strain partitioning structures, suggesting that they may share a common mechanism for their development (Figs. 3 and 5). Two segment boundaries, west Nepal–east Nepal and east Nepal–Bhutan, are located near the northward projection of Indian basement ridges (Gahalaut and Kundu, 2012) (Fig. 1), which may play a role in segmenting the range. The largest cluster of high-k_sn river channels is located in the east Nepal segment, which goes from the Thakkhola graben in the west to the Yadong-Gulu rift in the east, just over 500 km (Fig. 2). The high-k_sn cluster that defines the east Nepal segment is closer to the toe of the thrust wedge (60–70 km) than any other cluster overlying a weakly coupled portion of the MHT (Fig. 2). Only west Nepal and Bhutan contain high-k_sn clusters closer to the range front (30–50 km), but these overlie strongly coupled portions MHT (Fig. 2). In both west Nepal and Bhutan, a second high-k_sn cluster parallels the first but overlies the weakly coupled portion of the MHT 80–120 km from the toe of the thrust wedge. In west Nepal, this has been interpreted to reflect an active MHT flat connecting a hinterland ramp with the more foreland ramp (Harvey et al., 2015); our analysis suggests that a similar geometry may be present in Bhutan.

The segments documented here are typically 200–300 km in orogen-parallel length, with the exception of east Nepal which is 500 km long (Figs. 2 and 5; and Table 3). This regular spacing suggests that segment boundaries and Himalayan-Tibetan rifts may be the product of a frequency-dependent mechanism, possibly buckling of the downgoing Indian plate. Buckling of slabs is a function of arc curvature and slab thickness, which, in a homogenous medium, would result in a set wavelength and amplitude (Yin, 2000). Another possibility is that both rifts and k_sn segment boundaries are the geometric consequence of the curvature of the Himalayan arc (McCaffrey and Nabelek, 1998).

Structural analysis of the deformed Siwalik Group at the leading edge of the MHT in western Nepal indicates the presence of north-striking lateral ramps in the Main Frontal thrust (Mugnier et al., 1999) that may extend far into the thrust wedge as lateral ramps in the MHT. It is also observed that the intersection of Indian basement ridges and the Main Frontal thrust coincide spatially with the Uttarakhand–west Nepal and west Nepal–east Nepal segment boundaries, and project near the east Nepal–Bhutan boundary (Fig. 1) (Gahalaut and Kundu, 2012). The timing of plastic deformation and metamorphism varies across the Faizabad ridge (west Nepal–east Nepal boundary), suggesting that these differences could be explained by lateral ramps in the MHT produced by subduction of these features (Gibson et al., 2016). Based on our results, these hypothesized lateral ramps may also explain the margins of the strain accumulation segments. Another possibility is that lateral ramps are generated by slip on syncollisional normal faults striking at a high angle to the thrust wedge, such as the Thakkhola graben, Kaurik-Chango (Leo Pargil), Gurla Mandhata–Humla, and Yadong-Gulu rift systems (Fig. 1A). If true, slip on these extensional structures is coeval with thrusting below. Motion on these upper crustal structures and thrusting at lower crustal levels along lateral ramps along the MHT could produce orogen-normal shear zones that are kinematically linked with upper and lower structural levels, as has been previously suggested by Harvey et al. (2015).

Oblique convergence at subduction zones partitions strain into arc-normal and arc-parallel components (McCaffrey, 1992; Tikoff and Teyssier, 1994; Braun and Beaumont, 1995; Ellis et al., 1995; McCaffrey et al., 2000), and is present in the Himalaya as well. The most striking Himalayan example exists in the western Himalaya–Tibetan Plateau where the right-lateral Karakorum fault feeds slip through a major extensional stepover at the Gurla Mandhata core complex to the right-lateral Western Nepal fault system and then to the Main Boundary thrust near the toe of the thrust wedge (Armijo et al., 1986; McCaffrey and Nabalek, 1998; Murphy and Copeland, 2005; Styron et al., 2011; Murphy et al., 2014; Whipp et al., 2014; Silver et al., 2015). This linked system of faults divides the Himalaya into a western obliquely convergent sliver and a central orthogonally convergent region. At the eastern edge of the orthogonally convergent region is the Yadong-Gulu rift, which extends from central Tibet to the Lesser Himalaya and has been suggested to be a strain partitioning boundary (Wu et al., 1998; Velasco et al., 2007; Drukpa et al., 2006; Antolín et al., 2012), although the case for it is less clear than for the Karakorum fault–Western Nepal fault system. The Karakorum fault–Western Nepal fault system and Yadong-Gulu rift separate orthogonally convergent eastern Nepal from its obliquely convergent neighbors. Upon crossing either of these structures, the change is not simply kinematic; the spatiotemporal pattern of seismicity also changes. East Nepal features a continuous band of seismicity interpreted to mark a crustal scale ramp in the MHT (Pandey et al., 1995, 1999; Cattin and Avouac, 2000). Upon crossing from east to west the Western Nepal fault system, this band of seismicity splits, becoming two distinct bands (Cattin and Avouac, 2000; Harvey et al., 2015); similarly, upon crossing from west to east the Yadong-Gulu rift, the band of seismicity disappears, becoming a diffuse scattering of earthquakes across the Kingdom of Bhutan (Velasco et al., 2007; Drukpa et al., 2006) (Fig. 2). The strike-slip nature of the Western Nepal fault system separating west and east Nepal has been documented (Murphy et al., 2014; Silver et al., 2015); however, fault-slip data implicating strike-slip motion along the Yadong-Gulu rift is lacking. Nevertheless, the High Himalaya of Bhutan is offset from that of east Nepal and Sikkim by 70 km in a left-lateral sense (Wu et al., 1998); studies of seismicity find distributed mid-crustal strike-slip deformation across southern Bhutan (Velasco et al., 2007; Drukpa et al., 2006). Magnetic remanence in western Bhutan requires clockwise vertical-axis rotation, consistent with left-lateral slip on the Lingshi fault, an oblique-slip normal fault that runs parallel to the southern Yadong-Gulu rift on the Bhutanese side of the border (Antolín et al., 2012). The pattern of seismicity in Nepal is reflected in river channel steepness as a 500-km-long cluster of high-k_sn rivers located immediately north of the band of seismicity. In west Nepal there are two clusters of high-k_sn rivers, each associated with a band of seismicity. However, in Bhutan the pattern of river channel steepness is distinct from the diffuse distribution of earthquakes; high-k_sn rivers are clustered in two narrow bands, one near the Tibetan Plateau margin and one near the southern border with India. We interpret these two bands to reflect the presence of thrust ramp duplexes or antiformal stacks generated by tectonic wedging, upon which strain accumulates at a relatively slower rate than in east Nepal.

The High Himalaya overlies weakly coupled portions of the MHT (Stevens and Avouac, 2015) and separates a stably sliding décollement beneath southern Tibet from strongly coupled portions of the MHT, suggesting that the High Himalaya acts as a dashpot storing interseismic strain (Bilham et al., 1997, 2001; Jouanne et al., 2004; Feldl and Bilham, 2006; Meade, 2010; Ader et al., 2012). In this model, southernmost Tibet and the High Himalaya bend elastically to accommodate India-Asia convergence until the elastically stored energy is sufficient to overcome the static friction on the strongly coupled portion of the MHT. This buildup and release was highlighted by the 25 April 2015 Gorkha earthquake, which resulted in ∼0.5 m of surface lowering of the leading edge of the High Himalaya and 1 m of uplift of the Lesser Himalaya (Lindsey et al., 2015; Elliott et al., 2016). Despite its large magnitude (M_w 7.9), this earthquake did not rupture the strongly coupled MHT. Rather, it propagated along the boundary between weakly and strongly coupled areas, uplifting the intervening moderately coupled region and transferring strain equivalent to a M_w 6.9 earthquake from the High Himalaya to the frontal portion of the range (Grandin et al., 2015).

The spatial distribution of historical earthquake ruptures in the Himalaya is an area of ongoing research (Ambraseys and Jackson, 2003; Kumar et al., 2006, 2010; Mugnier et al., 2013; Murphy et al., 2014; Berthet et al., 2014; Bollinger et al., 2014; Rajendran et al., 2015; Hossler et al., 2016). Our analysis reveals a rough correspondence between ruptures and segment boundaries delineated by high-k_sn clusters. While major earthquake ruptures are each largely confined to a single segment, the edges of rupture patches in some cases bleed over into adjacent segments. This indicates that segmented strain accumulation controls where and when major earthquakes initiate, but that segment boundaries are not necessarily barriers to rupture propagation. Rupture patches of large Himalayan earthquakes since A.D. 1900 correspond to segments defined by high-k_sn clusters. In our study area, there have been seven major Himalayan earthquakes during the past 500 yr, of which six are confined to individual segments. Only the A.D. 1505 M_w 8.2 earthquake is an exception; it ruptured the entire west Nepal segment, the eastern edge of the Uttarakhand segment, and the western edge of the east Nepal segment.

The most striking correspondence between the k_sn analysis and the spatiotemporal distribution of historical earthquakes lies within the Nepal-Bhutan segments. The west Nepal segment has not ruptured since 1505 (Ambraseys and Jackson, 2003) and Bhutan has not ruptured since 1713 (Berthet et al., 2014), while the east Nepal segment has ruptured three times over the last three hundred years in 1833, 1934, and 2015 (Fig. 6). This apparent dichotomy may simply reflect a lack of resolution in our knowledge of the Himalayan seismic cycle, or it may indicate that the central Nepal segment accumulates and releases strain more efficiently than its neighbors. This latter idea is supported by geomorphologic and kinematic considerations.

The east Nepal strain accumulation segment occupies the most orthogonally convergent sector of the orogen (Figs. 3 and 5), consistent with pure shear shortening, and hosts a cluster of high-k_sn channels that is broader, longer, and closer to the range front than those anywhere else in the Himalaya. The east Nepal segment contains eight out of the ten Himalayan peaks higher than 8000 m, and has well developed frontal anticlines that collect all of the High Himalayan catchments into three outlets along the range front. Taking the highest topography in the range as a measure of accumulated High Himalayan strain, and the most well-developed frontal anticlines in the range as a measure of Sub-Himalayan strain accumulation, suggests that range-normal strain accumulation and release is more efficient here than anywhere else in the Himalaya.

West Nepal and Bhutan lie within the obliquely convergent western and eastern sectors of the Himalaya, respectively, in which some portion of India-Asia convergence is partitioned into orogen-parallel motion (Drukpa et al., 2006; Velasco et al., 2007; Styron et al., 2011; Antolín et al., 2012; Murphy et al., 2014; Whipp et al., 2014). Changes in the size of the hinterland region accumulating strain, its proximity to the toe of the thrust wedge, and lateral strain partitioning provide a plausible explanation for the seemingly disparate seismic cycles of west Nepal, east Nepal, and Bhutan.

Looking at earthquake ruptures even further back, we encounter an apparent change in seismic behavior. There have been four surface-rupturing events identified between A.D. 1000 and A.D. 1500: 1100, 1255, 1344, and 1430, three of which have inferred rupture lengths between 500 and 1000 km requiring magnitudes of M_w 8–9 (Sapkota et al., 2013; Berthet et al., 2014; Mugnier et al., 2013). While this may be the result of conflating two (or more) smaller ruptures spaced closely in time, it may be evidence of a >1000 yr seismic cycle that culminates in an ∼M_w 9.0 earthquake (Kumar et al., 2006; Stevens and Avouac, 2016). Viewing the recent sequence of large earthquakes in central Nepal (in 1833, 1934, and 2015) in this context, only the 1934 surface-rupturing event represents a significant decrease in accumulated strain. The 1833 and 2015 earthquakes instead represent strain transfer from the hinterland to the foreland. This idea is supported by an apparent strain deficit between the magnitude of convergence during the past 500 yr and the magnitude of strain released in great earthquakes, which is estimated to be equivalent to four M_w 8 earthquakes or a single M_w 9 earthquake (Bilham and Ambraseys, 2005; Meade, 2010; Stevens and Avouac, 2016).

The drainage network of major Himalayan rivers is segmented similarly to channel steepness (Fig. 1), suggesting that channel steepness and drainage network patterns form in response to the same tectonic forcings. The Nepal segments are characterized by large catchments in which numerous parallel drainages are collected into a few pour points at the range front. This process is driven by the location of active anticlines in the Sub-Himalaya, which grow more rapidly than the individual rivers can incise through them, producing lower river reaches that flow parallel to the range front for 70–100 km. The result is a 760-km-long segment of the range hosting only three pour points entering the foreland basin—the Karnali, Narayani, and Arun rivers. In west Nepal, near the Tibetan Plateau boundary, in the upper reaches of the Karnali river, the Mugu Karnali in the west and Dolpo Karnali in the east flow along the axis of the Dolpo syncline, paralleling the range front for 100 km and 150 km, respectively. Axial drainage along the leading edge of the growing Dolpo anticline (Cannon and Murphy, 2014) by the Karnali’s two main trunk channels is funneled into a zig-zag drainage pattern in one of the largest erosional bites in the Himalaya (Fig. 1). The lower reaches of the Karnali contain two hairpin turns, suggesting recent drainage network reorganization by the development of active anticlines at the range front. The 270-km-long Bhutan segment contains five catchment basins, the largest of which hosts a parallel drainage network characteristic of steep, active mountain fronts. Tectonic activity along the range front in Bhutan is expressed differently than in the western and central Himalaya. To the west, the range front in Nepal is characterized by closely spaced anticlines separated from the range to the north by broad flat valleys (duns). This is in stark contrast to Bhutan, where foothill anticlines are intermittent, duns are narrow, the Main Frontal thrust is difficult to recognize, and the first topographic step is commonly associated with the Main Boundary thrust. However, these differences are minor compared to the presence of the Shillong Plateau in the eastern Himalaya foreland basin. The Shillong Plateau is located 100 km south of the range front and has >1 km of relief. It has been characterized as a popup structure bounded to the north and south by oppositely dipping thrusts, however its significance for the thrust wedge and its implications for the interaction between climate and tectonics are debated (Adams et al., 2015; Bilham and England, 2001; Clark and Bilham, 2008; Grujic et al., 2006; Najman et al., 2016). Whatever role the Shillong Plateau plays in accommodating India-Asia convergence, it seems clear that at the range front strain is accumulating at a lower rate than in Nepal, as all drainages flow across the trace of the Main Frontal thrust without disruption.

Despite numerous geologic and geophysical studies of the Himalayan plate boundary, basic questions remain unanswered. Foremost is, what is the geometry of the subduction interface? Seismic reflection surveys have resolved ramps and flats in the MHT (Hauck et al., 1998; Caldwell et al., 2013), although seismic refraction surveys have failed to do so (Nábělek et al., 2009). Geodesists are able to model the pattern of MHT fault coupling using a single inclined plane by assuming that all significant variability in the position and dip of MHT ramps is located within strongly coupled portions (Bollinger et al., 2004b; Jouanne et al., 2004; Stevens and Avouac, 2015). Nevertheless, there is evidence for the existence of MHT ramps located within the weakly coupled portions of the MHT. Balanced cross sections (Schelling and Arita, 1991; DeCelles et al., 1998; Robinson, 2008; McQuarrie and Ehlers, 2015; Hubbard et al., 2016), electrical conductivity (Lemonnier et al., 1999), tectonic geomorphology (Seeber and Gornitz, 1983; Cannon and Murphy, 2014; Harvey et al., 2015; Morell et al., 2015), microseismicity (Fig. 2) (Pandey et al., 1999), and Bayesian modeling of possible MHT geometries (Elliott et al., 2016) all suggest the presence of thrust ramps in essentially the same locations, at the boundary between the strongly and weakly coupled portions of the MHT. Changes in the position of high-k_sn clusters are most pronounced between west Nepal and east Nepal, where high-k_sn clusters overlying weakly coupled portions of the MHT are offset by 60 km along a sharp boundary coincident with the Gurla Mandhata–Humla core complex (Harvey et al., 2015) (Figs. 2 and 5). This indicates that the geometry of the MHT changes along strike. That the largest offset of high-k_sn clusters occurs here has kinematic significance. The Gurla Mandhata–Humla core complex is thought to be an extensional stepover feeding slip from the right-lateral Karakorum fault in southern Tibet into the right-lateral cross-orogen West Nepal fault system (Murphy et al., 2014; Silver et al., 2015).

East Nepal is located in the most orthogonally convergent part of the Himalayan arc (Fig. 5) and during that past 500 yr has had the most rapid recurrence of M_w >7 earthquakes anywhere in the range (Fig. 6). Convergence obliquity and great earthquake recurrence times increase to the west and east. How much convergence obliquity modifies the geometry of the plate margin or modulates its seismic behavior remains unknown, but is increasingly seen as an important seismotectonic variable (McCaffrey and Nabalek, 1998; Jouanne et al., 1999; Banerjee and Bürgmann, 2002; Styron et al., 2011; Murphy et al., 2014; Whipp et al., 2014; Kundu et al., 2014).

Middle Miocene to recent internal deformation of the Himalayan thrust wedge appears to be dominated by thrust duplexing (Bollinger et al., 2004a; Herman et al., 2010; Webb, 2013; He et al., 2015; Larson et al., 2015; Yu et al., 2015). At the surface, duplexing is characterized by a high ratio of crustal thickening to surficial horizontal shortening. This criterion is met in the hinterland of west Nepal where the broad Dolpo anticline records vertical thickening to surficial horizontal shortening ratios of between 2:1 and 5:1, depending on the thickness of the repeated layer (Cannon and Murphy, 2014). Duplexes typically form by underplating, whereby slices of the footwall are accreted to the hanging wall (Fig. 7), however no structural window in the High Himalaya is documented that provides an unambiguous view of the mechanism thickening the hinterland. Alternatively, in the plastically deforming lower crust, repeated tectonic wedging (Webb et al., 2007; Yin, 2006) at the BPT could produce a duplex restricted to the hanging wall of the MHT (Fig. 7). Tectonic wedging focused at the BPT could alleviate some of the slip deficit by allowing the lower crust to accommodate more convergence than the upper crust. In this view, hinterland thickening results from repeated insertion of plastic wedges at a ramp where the MHT BPT is localized, forming an antiformal stack. As the active wedge is inserted at the base of the stack, it forces the preceding wedges up, building topography and steepening river channels at the surface. The brittle upper crust overlying the antiformal stack absorbs elastic strain imparted by the south-southwest–directed sliding of southern Tibet recorded by GPS (Ader et al., 2012). Hinterland lower to mid-crustal duplexing provides a plausible explanation for the coincidence of interseismic strain accumulation in the thrust wedge, spatial clusters of high k_sn values, and the highest topography in the range without calling upon a delayed erosional response to ramp migration (Grandin et al., 2012), and diminishes the need for an unbalanced earthquake cycle (Bilham and Ambraseys, 2005; Meade, 2010; Stevens and Avouac, 2016) by allowing the lower crust to accommodate more permanent strain than the upper crust. Recently, several mid-crustal discontinuities have been identified within the Greater Himalaya sequence of Nepal predating the main period of the linked Main Central thrust–South Tibetan detachment motion in the late early to middle Miocene (Carosi et al., 2010; Larson et al., 2015; Larson and Cottle, 2014; Montomoli et al., 2013). These discontinuities have been interpreted as shear zones bounding horses within a hinterland duplex that was exhumed during thrusting on the Main Central thrust. Hinterland duplexes create permanent uplift by thickening the lower crust while storing elastic strain in the upper crust, the net result of which is collocated increased rates of rock uplift and elastic strain accumulation. Currently, we cannot distinguish between these two duplexing models, although the presence of isotopic signatures in anatectic melts of the Gurla Mandhata core complex provides compelling evidence for migration of Indian-affinity material across the MHT (Murphy, 2007).

Rapid rates of inferred rock uplift in the Himalaya, recorded by k_sn analysis, are predominantly localized within the High Himalaya above weakly coupled portions of the MHT. The time scales required to encode k_sn signals in the landscape (10⁵ yr) provide an intermediary time step between GPS data sets and low-temperature thermochronology. Clusters of high-k_sn rivers are offset from each other in map view across Himalayan-Tibetan strain partitioning structures and rift systems. This suggests that the location of linked strike-slip and extensional structures is fundamental to the structural segmentation of the orogen. Regions of rapid rock uplift in the High Himalaya separate stably sliding southern Tibet from the Himalayan seismogenic zone, and acts as a dashpot accumulating and storing elastic strain. In addition to acting as a reservoir for accumulating elastic strain, collocation of clusters of high-k_sn streams with the highest topography in the Himalaya indicates that they also record crustal thickening. As such, clusters of high-k_sn streams represent an accumulation of combined interseismic and coseismic strain. Structural segmentation of the Himalaya based on spatial clustering of high k_sn values reveals an association with the spatiotemporal distribution of historical great earthquakes in the central Himalaya. This relationship suggests that the density and spatial dimensions of regions characterized by continuous high values of k_sn, their distance from the toe of the thrust wedge, and the degree of convergence obliquity impact the recurrence of major earthquakes in active orogenic systems.

The manuscript benefited from reviews by Dr. Jon Harvey and an anonymous reviewer. The authors would especially like to thank Dr. Harvey for thorough and insightful reviews which led to substantiative improvements to the manuscript.