Recent uplift of Chomolungma enhanced by river drainage piracy

Chomolungma, also known as Mount Everest or Sagarmāthā (Fig. 1a), is the highest mountain on Earth, towering over other peaks in the Himalayas (Fig. 1b). At 8,849 m above sea level, Chomolungma is ~250 m higher than the other tallest peaks in the Himalaya. This is surprising given the relative along-strike uniformity of tectonics in the Himalayas, which provides mountain peak buoyancy with local, small-scale variability^1,2, and relatively uniform climatic conditions and erosional processes^3,4. Despite this, Chomolungma deviates from the linear trend of the elevation of Himalayan peaks plotted against their elevation rank (Fig. 1b). A simplistic analysis hints that this anomalous elevation (on the order of one to two hundred metres) could be the result of a scenario where peak elevations were relatively uniformly distributed, but an additional component of elevation has been added to only some of the peaks (including the highest peak). Only the highest of these affected peaks would disrupt the linear trend, whereas lower peaks affected by the additional elevation simply move up the ranking (Fig. 1b). Further evidence suggests a disconnect between the short-term rock uplift rate of Chomolungma measured from GPS data (~2 mm yr⁻¹)⁵ and its long-term rock uplift rate recorded by thermochronology (~1 mm yr⁻¹)⁶ (Supplementary Text 1), which raises the question, is there an underlying mechanism raising Chomolungma’s anomalous elevation ever higher?

a, A photograph of Chomolungma taken from its northern side. At this ___location, the Arun River flows towards the north away from Chomolungma before turning south through the Arun Gorge. b, Mountain peaks above 7,500 m in the Himalaya ranked by elevation. Chomolungma is off the almost linear relationship between rank and elevation. This relationship can be approximated with peaks drawn from a uniform distribution of elevations between 7,500 m and 8,400 m (blue curve) with an additional 300 m of uplift affecting 30% of these peaks, including the highest point (red curve). c, The enigmatic drainage network of the Kosi River. Local relief is highest where the gorge cuts through the Himalaya and elevation varies from 8,481 m to 1,140 m over a distance of only 35 km. The main faults of the Himalaya are the Main Frontal Thrust (MFT), the Main Boundary Thrust (MBT), the Main Central Thrust (MCT) and the South Tibetan Detachment System (STDS). The zone of high rock uplift rates is associated with the development of a mid-crustal duplex and is resolved with thermochronometry⁶. The rock uplift rate in this zone is set to 4 mm yr⁻¹ and is 0.5 mm yr⁻¹ to the south and 1 mm yr⁻¹ to the north. d, The χ map for the Kosi River catchment. e, Most of the divides between the catchments of tributaries and across the Himalayan divide are relatively stable. f, In contrast, in the Arun River and the tributaries to the east and west, χ varies considerably.

Full size image

One potential mechanism to explain Chomolungma’s anomalous elevation and recent increased rock uplift rates becomes apparent when inspecting its surrounding rivers. The Arun River, a major tributary of the Kosi River, drains a large area to the north of Chomolungma before turning south, passing by the world’s tallest peak and cutting a deep gorge through the core of the Himalayas. Ever since the earliest expeditions to Chomolungma, there has been debate as to whether this unusual drainage pattern represents antecedent drainage patterns that the Himalayas have pushed through⁷ or if the Arun River has captured Tibetan rivers⁸. This capture event could be related to localized tectonics triggered by fluvial erosion^9,10, though the timing of this event is unconstrained.

This scenario of recent drainage piracy has implications for Chomolungma’s elevation. Drainage piracy, facilitated by headward erosion, alters water flow and sediment transport routes by reshaping the drainage divide and network^11,12,13,14. If recent drainage piracy is responsible for the drainage pattern of the Kosi River and Arun River, an increase in river discharge due to capture would drive gorge incision. This incision would trigger isostatic compensation, causing Earth’s crust to rebound, resulting in surface uplift of the unincised parts of the surrounding area, including the mountain peaks^15,16. The size of the surrounding area that rebounds due to incision is controlled by the strength of the lithosphere. Flexural isostatic rebound, and resulting surface uplift, in response to capture-driven incision of the Arun River could therefore explain part of the observed elevation anomaly of Chomolungma.

The identification of drainage piracy mainly relies on morphological and geological clues. Drainage piracy can be inferred from large variations in χ values (a parameter that characterizes river network topology and geometry) on either side of the divide¹⁷, or the presence of wind gaps within divides^18,19, but neither indicator provides absolute time constraints. In contrast, provenance changes of sedimentary basins provide indirect insights on drainage reorganization and potential timing constraints^20,21. Direct time constraints can be inferred by assessing the rapid incision of rivers²², but isolating the signal of capture from changes in tectonics or climate is challenging.

To understand the drainage evolution of the Kosi River, and the implications of this evolution for Chomolungma, we combine a simple stream power model (SPM) with the Neighbourhood Algorithm²³. This approach allows us to reconstruct the capture process of a river, estimate fluvial erosion parameters and quantify uncertainties. Our study sheds light on recent rock uplift rates and the anomalous elevation of Chomolungma, unveiling the interplay between river capture and mountain uplift.

Located in the heart of the central Himalayas (Fig. 1c), the Kosi River, and its major tributary the Arun River, traverse the highest regions of the Himalayas. The Arun River drains the southern expanse of Tibet and the northern slopes of Chomolungma before flowing through a narrow gorge that drops 7 km in elevation over 35 km (Fig. 1c and Extended Data Fig. 1).

The Kosi River has a distinctive river profile and planform characteristics (Fig. 2a and Extended Data Fig. 2). Whereas there are other rivers in the Himalayan region that have east–west headwaters, such as the Sutlej and Ghagra (Extended Data Fig. 2a), their longitudinal profiles are typically concave. In contrast, the longitudinal profile of the Kosi River shows a pronounced convexity (Extended Data Fig. 2h). The modern topography provides compelling evidence that the Kosi River system is in an unsteady state in which erosion rates are not balanced by rock uplift rate (Fig. 1d–f).

a,b, Distance to base level versus elevation (a) and χ versus elevation (b) for the Arun River and rivers to the east (Tamor River) and west (Dudh Kosi River and Likhu Khola River). The overall form of the rivers is similar in a but the Arun River has the additional area from the Tibetan Plateau. The transformation from distance to χ should collapse rivers that have experienced a common rock uplift rate on to the same curve. This is the case for the Tamor River, Dudh Kosi River and Likhu Khola River, and we argue these are in steady state with the regional rock uplift rate pattern. However, the central trunk channel does not collapse on the same overall curve.

Full size image

Examining a χ–elevation plot^24,25 of the region (Fig. 2b), it is clear that the eastern and western tributaries of the Kosi River have a similar rock uplift history. They converge on a common form and their tributaries exhibit similar behaviour. In contrast, the Arun River, the major tributary in the middle, deviates from this form, despite the similar lithology and tectonic setting. Whereas there is some evidence for differential tectonics, with thrust faults on either side of the Arun River, these tectonic processes are minor. We therefore suggest that the χ–elevation relationships of the Arun River are different due to an ongoing adjustment of river elevation in response to a recent increase in stream power (Extended Data Fig. 3).

Our approach focuses on modelling the capture process through changes in the upstream reaches. Before the capture event, we assume that all rivers are at topographic steady state (Supplementary Text 1) with the long-term rock uplift rate inferred previously with thermochronometry⁶, and the Arun River does not extend upstream of the capture point (Extended Data Fig. 4a–c). After capture, the upstream area of the modern Arun River is included, assuming it has its current morphology (Extended Data Fig. 4d). The tributaries to the east and west of Arun River do not experience this same increase in upstream drainage area and are assumed to remain in a state of equilibrium. In contrast, the Arun River responds to the addition of new upstream drainage area by capture, the timing of which is unknown and is referred to as t_c.

We simulate this capture event using the detachment-limited SPM²⁶:

where h is elevation, t is time and thermochronological data and thermokinematic modelling⁶ are used to constrain the rock uplift rate (U) (Extended Data Fig. 5). The upstream drainage area (A) evolves from a pre-capture network to a post-capture network at t_c, and S is the channel slope. The unknown model parameters are the area exponent (m), the slope exponent (n), the erodibility (K) and the time of capture (t_c). We model the Kosi River network, including rivers that are impacted by the capture event and stable rivers not impacted by the capture event. In our forward model, we use an implicit method^27,28 to predict river nodes through time. The initial condition assumes a topographic steady state landscape that is predicted using defined values of m, n and K. Using common values of n = 1, m = 0.45 and K = 8 × 10⁻⁶yr⁻¹, about 76 thousand years is required to erode the Arun Gorge from the initial topographic steady state landscape.

However, SPM parameters are uncertain, with some studies based on cosmogenic nuclide datasets suggesting that the slope exponent n > 2 (refs. ^29,30). Therefore, to propagate uncertainties from these parameters to uncertainties in the capture time, we use an inverse model. The Neighbourhood Algorithm²³ is applied to infer values of m, n, K and t_c. At each iteration of this algorithm, the selected values of m, n and K are used to build an initial, pre-capture landscape, that is, considered as an initial condition. The capture event is then simulated with these parameters and the capture time. Misfit is calculated for both the topographic steady state landscape not impacted by the capture event and the Arun River that is adjusting to the capture event. The ability of the model and data to resolve these model parameters is discussed in the Methods. A similar inversion scheme has been introduced by ref. ³¹ that focused only on fitting a single knickpoint, however, our approach includes simulation of the entire drainage network, providing improved resolution.

Modelling results highlight excellent agreement between model predictions and topography. Residuals are centred around zero (Extended Data Fig. 6), although some areas show notable differences between observed and predicted elevations. In general, the trunk streams are accurately predicted, whereas certain headwaters of tributaries may be susceptible to over- or under-prediction (Fig. 3a). We observe that rivers tend to be over- and under-predicted across drainage divides³² where active migration is evident, a factor not accounted for in our model. Furthermore, the highest parts of our river network are also poorly reproduced, which may be due to glacial erosion. Nonetheless, our model and its associated parameters show a good fit to the data, providing confidence in our results.

a, The model residuals show relatively low misfits across all aspects of the model ___domain. In some upper tributaries, residuals are large, and we attribute this to glacial erosion and ongoing drainage divide migration. b, The recovered profile matches the trunk stream well, and the overall form of the tributaries are reproduced. Furthermore, concavities on the tributaries are reproduced, highlighting the transient wave of incision associated with increased upstream drainage area. The zone of high rock uplift rate for the main stem of the Kosi River is shown in grey shading. c, The four model parameters are well resolved as shown by the tight clustering of low misfit values. The parameters m and n covary with one another, as would be expected; however, a clear minimum corresponding to m = 0.43 and n = 0.97 is observed. The star shows the ___location of the optimal solution. d, The parameters K and t_c also covary, but a best fitting value of K = 1.16 × 10⁻⁵ m^1–2m yr⁻¹ and t_c = 89 ka can be identified. The star shows the ___location of the optimal solution. We note that the units of K are variable because m varies. For the best fitting model parameter combinations of n = 1 and m = 0.42, the units correspond to m^0.16 yr⁻¹.

Full size image

Our best-fit values of m, n and K are 0.43 and 0.97 and 1.16 × 10⁻⁵ m^1–2m yr⁻¹, respectively (Fig. 3b,c), and the posterior probabilities are 0.45 ± 0.083, 1.07 ± 0.17 and −5.07 ± 0.58 for m, n and log₁₀ K, respectively. The posterior probability distributions for these parameters are shown in Extended Data Fig. 7. Although the ratio between m and n is almost identical to the value constrained for the Himalayas in previous studies^33,34,35,36, our value for n supports evidence from studies of fluvial erosion that river incision rate scales linearly with slope (n = 1)^37,38,39. This contradicts a number of studies suggesting that n may be greater than 1 in tectonically active regions^29,40 and in the Arun River⁴¹. However, these studies are based on catchment-wide erosion rates from cosmogenic nuclide concentration measurements that are sensitive to both fluvial and hillslope erosion and are based on a steady state assumption. In the Himalayas, for example, the transient wave of incision we identify would invalidate this key assumption. Our results also help explain some of the very high rates of erosion recorded in the catchment-wide cosmogenic nuclide data⁴¹. In particular, erosion rates up to 10 km per million years were attributed to local landsliding and deemed unreliable; however, such high transient rates are actually predicted by our model.

Our best-fit model estimates that capture of the Arun River occurred 89 thousand years ago (ka), whereas other well-fitting models suggest capture events between 50 ka and 100 ka (Fig. 3d). The posterior probability of the capture time is asymmetric with a maximum value of ~70 ka (Extended Data Fig. 7). This closely matches the age of pluvial sediments upstream of the capture point, dating back to about 30–90 ka (refs. ^42,43). It also supports the validity of our inferred capture timing, indicating that the pluvial sediments represent the final stages of a High Himalaya drainage divide rather than a recent blockage event caused by phenomena such as landslides.

In other Himalayan rivers, it has been suggested that catastrophic erosion events are essential to maintain high erosion rates^44,45 and cut deep gorges that are in topographic steady state with rock uplift rates. However, the similarity in our capture time to pluvial lake sediments and the landscape’s adjustment to increased drainage area suggests additional catastrophic events (failure of small blockages of the narrow gorge due to landslides or glaciers) might not be necessary for the Arun Gorge formation. Furthermore, the Arun River’s morphology supports ongoing landscape adaptation to the larger upstream drainage area (Figs. 1d and 2).

It is important to recognize that we do not account for the spatial variability of precipitation. Modern precipitation patterns across the Himalayas show an increase from south to north before decreasing within the southern Tibetan Plateau. Taking this into account would imply that capturing the Tibetan portion of the Kosi catchment, where precipitation rates are lower, would result in a smaller change in discharge along the Kosi River than simulated in our model. As a result, the time required to produce the observed changes in headwater stream elevations and the retreat of knickpoints along its tributaries may be even longer than our estimate. It is worth noting, however, that there is limited understanding of how precipitation evolved during the Pleistocene, how the landscape changes we have identified would affect precipitation patterns, and the precise impact of precipitation on erosion rates. Therefore, we prefer our simple model, in which the capture time corresponds to the age of pluvial sediments, coupled with qualitative analysis of extreme events and the influence of precipitation.

Following the capture event resolved by our model, downstream incision rates increase due to the increased upstream drainage area. Here we define incision rate as the difference between the rock uplift rate and the river erosion rate. In this way, incision is a non-steady process that generates topographic relief. According to our model results, downstream of the capture point, incision rates vary between 5 and 10 mm yr⁻¹ in the upper reaches, 1–5 mm yr⁻¹ in the middle reaches and less than 1 mm yr⁻¹ in the lower reaches (Fig. 4a). This result is consistent with previous incision rate estimates based on terrace dating near Tumlingtar⁴⁶.

a, The difference between rock uplift rate and erosion rate is greatest in the High Himalaya. Due to the recent change in fluvial incision, it is expected that the highest parts of the landscape have not responded to this increased erosion. Therefore, these areas will experience isostatic surface uplift. b, The flexural response to enhanced erosion of the Kosi River leads to enhanced surface uplift across the highest Himalaya for T_e = 20 km.

Full size image

The incision of the Arun Gorge (Extended Data Fig. 8), driven by the capture of the upstream drainage area, represents a dynamic and non-steady component of erosion. In the main rivers flanking the Arun River to the east and west, the high rates of erosion balance the rock uplift and peaks remain at a constant elevation (Fig. 4a). Conversely, the transient component of incision in the Arun River drainage triggers an isostatic response. Due to flexure of the lithosphere, an area extending over tens of kilometres experiences renewed surface uplift, contributing to the elevation of the world’s highest mountains. Montgomery¹⁶ suggested that 20–30% of Himalayan peak elevation can be attributed to isostatic rebound induced by river incision. This conclusion extends the earlier findings of Wager¹⁵. However, these results were unusually high because they assumed that an extensive flat topography existed before all the material between the peaks was eroded. On the basis of normal response times of landscapes⁴⁷, and given that the mean elevation has been consistent in the Chomolungma area since the Miocene⁴⁸, we expect this area has achieved topographic steady state before Quaternary fluvial capture. Following capture, enhanced incision produced an isostatic component of rock uplift (Extended Data Figs. 9 and 10). On the basis of our inferred SPM parameters, the wave of fluvial incision associated with the river capture has not propagated significantly along the tributaries. Under topographic steady state assumptions, the surrounding peaks are therefore expected to be eroding at the previous rock uplift rate and so experience flexural isostatic surface uplift. With an effective elastic thickness of 10–30 km, this phenomenon would lead to an increased surface uplift rate for Chomolungma of 0.16–0.53 mm yr⁻¹. This isostatic response contributes an additional rock uplift component equivalent to approximately 10–50% of the current GPS-derived rock uplift rate⁵. GPS stations in proximity to the area of rapid non-steady incision will be more strongly affected, with important implications for the analysis of interseismic strain accumulation and earthquake hazard⁴⁹.

We suggest that Chomolungma partly owes its anomalous elevation to drainage capture and river incision (Fig. 4b), though tectonic processes governing crustal thickness and rock uplift are clearly the fundamental cause of Chomolungma’s elevation⁴. Our study uncovers a previously unrecognized, additional mechanism of rock uplift active since river capture. River incision lowers the elevation along the river channel, but the regional area is subject to rock uplift by isostatic rebound. Because tributary rivers take time to respond to this incision, erosion rates lag behind the isostatic rock uplift rate leading to surface uplift. To estimate the magnitude of additional elevation, we calculate the total amount of river incision between the capture time and the present. This provides an estimate of the eroded volume, which can be converted to flexural isostatic surface uplift. In particular, the integrated isostatic surface uplift of ~15–50 m may help explain part of the anomalous height of Chomolungma (8,849 m) with respect to K2 (8,611 m) and the next three highest peaks in the Himalaya (8,586 m; 8,516 m; 8,485 m), which differ in elevation by no more than 100 m between one peak and the next. However, Chomolungma and the neighbouring highest peaks (Lhotse and Makalu) experienced similar amounts of surface uplift, moving up the rankings, but not disrupting the linear relationship between peak elevation and ranking (Fig. 1b). In contrast, Chomolungma has been raised to an even higher elevation. We therefore attribute part of Chomolungma’s anomalous elevation to the enhanced surface uplift it experiences due to Arun River gorge incision.

Forward and inverse models

We applied a series of forward models to simulate river evolution, using the stream power model to predict river network elevations based on model parameters and rock uplift rates. We assumed a temporally homogeneous but spatially variable uplift rate. The rock uplift rate is based on the extensive thermochronometric studies carried out in this area, which has constrained a narrow band of enhanced rock uplift rate associated with movement over a deeper ramp structure⁶ (Extended Data Fig. 5). In the initial condition of the model, the Arun River contains only the portion downstream the capture point, whereas the river morphology of other tributaries (Tamor River, Dudh Kosi River and Likhu Khola River) is consistent with the current state (Extended Data Fig. 4c). The capture point is determined by inspection of the digital elevation model (DEM), but our results are robust with respect to the exact capture point. Considering that these tributaries do not have obvious knickpoints and the channel profiles are essentially in equilibrium, we have assumed that the rivers are in topographic steady state at the start of the model (Extended Data Fig. 4a,b). After capture, the morphology of the western and eastern tributaries remains unchanged and therefore at topographic steady state. The Arun River adds an upstream portion that is consistent with the present morphology (Extended Data Fig. 4d) and begins to incise due to the increase in upstream area.

We used TopoToolbox 2⁵⁰ to extract the drainage network from NASADEM elevation data⁵¹ at a resolution of ~ 30 m and resampled to 100 m to improve computational speed. We defined the river network as all areas draining more than 1 km². The residuals in Fig. 3a and Extended Data Fig. 6 are calculated based on the difference between the model and real elevations at the river nodes.

Considering that these parameters have been reported to be highly variable in the natural environment, we need to propagate the uncertainties on m, n and K to capture time. The nonlinear inverse model uses the neighbourhood algorithm (NA)²³ to propagate uncertainties and search optimal model parameters. We solve for four free parameters—the area exponent (m), the slope exponent (n), the erodibility (K) and the time of capture (t_c). The m value is between 0.2 and 0.7, n from 0.5 to 2, K from 1 × 10⁻⁷ to 1 × 10⁻⁴ m^1–2m yr⁻¹ and t_c from 0 to 500 ka.

It is important to consider the potential of the model and data to resolve the model parameters. The model and data are able to resolve the area exponent (m) from the elevation and fluvial response times of river nodes along stable river profiles and the requirement that these nodes converge along a consistent elevation–response time curve⁵². The value of slope exponent (n) is resolved using a similar principle, taking into account the shapes of tributaries to the Kosi River, which respond to the same wave of incision but may have different local gradients. Erodibility (K) is determined by ensuring that the stable rivers match the thermochronologically constrained rock uplift rate. It is important to note that m, n and K are interrelated⁵³. However, by incorporating transient river profiles, knowledge of spatially variable rock uplift rates and a relatively simple geomorphic model, we can resolve these parameters independently. The relative timing of capture is determined by assessing the differences between the stable pre-capture Kosi River and its tributaries compared with the modern transient profile. The ability to resolve K impacts our ability to determine the absolute value of t_c; however, we find that these parameters are independently resolved.

The algorithm works by running a series of forward models that predict observed data. We use an L1 norm to calculate a misfit (M) between the model predicted elevations and the observed elevations.

where N is the number of river nodes, p_i is the predicted river node elevation and o_i is the observed river node elevation. This ensures that our parameter estimation is robust with respect to outliers. Only river nodes of Tamor River, Dudh Kosi River, Likhu Khola River and downstream of the capture point for Arun catchment are used to assess the performance of the forward model and model parameters. River nodes with elevations higher than 4,000 m were excluded to avoid the expected impact of glacial erosion; river nodes with upstream drainage areas smaller than 5 × 10⁶ m² were excluded to ensure that we do not include hillslopes. A total of 51,000 forward models were run split into 51 iterations with 1,000 models per iteration. At each iteration, the NA resamples promising areas of parameter space but maintains flexibility to explore and discover new promising areas.

To provide robust parameter estimation, we integrate information across parameter space using the second stage of the NA. The posterior probability density functions are based on the likelihood function L, which is defined as:

where σ is the data uncertainty (10 m) and M is misfit. The 1D marginal of the posterior probability distributions for the four parameters are shown in the Extended Data Fig. 7.

Flexural model

After the river capture, the isostatic rebound induced by non-steady incision would produce an additional surface uplift. Here we used a two-dimensional flexure model⁵⁴ to calculate the vertical deflection of an elastic by a Fourier transform approach. In this model, the crust is considered as an elastic sheet with a flexural rigidity (D).

where E is Young’s modulus (1 × 10¹¹N m⁻²), v is Poisson’s ratio (0.25) and T_e is effective elastic thickness.

On the basis of the best-fit model, we calculated present-day isostatic uplift rates (Fig. 4b and Extended Data Fig. 9) according to the present-day incision rate (Fig. 4a). We also estimated the total amount of isostatic uplift since the timing of capture (Extended Data Fig. 10) based on the total non-steady incision amount (Extended Data Fig. 8). On the basis of seismological data and gravity data, the effective elastic thickness of Himalaya is 10–30 km (refs. ^55,56) (Supplementary Text 2). Our analysis explored a reasonable range of the effective elastic thickness (T_e) at 10, 20 and 30 km. The maximum rate of current surface isostatic uplift induced by river incision ranged from 0.2 to 1.1 mm yr⁻¹ (Extended Data Fig. 9). Specifically, Chomolungma experiences uplift rates of 0.16 to 0.53 mm yr⁻¹, whereas Makalu exhibits uplift rates of 0.18 to 0.77 mm yr⁻¹. Since the capture, the maximum total uplift due to non-steady incision ranges from 16 to 84 m (Extended Data Fig. 10). This includes 13 to 41 m uplift for Chomolungma and 14 to 60 m for Makalu. We note that we do not model hillslope erosion. We expect that hillslopes immediately adjacent to the incising channel will also experience enhanced erosion. This would increase the area of enhanced erosion and increase the isostatic response. In this way, our isostatic calculations are minimum estimates of surface uplift rates and total surface uplift.