Formation of Longitudinal River Valleys and the Fixing of Drainage Divides in Response to Exhumation of Crystalline Basement

Variations in rock strength act as a first‐order control on mountain landscapes. However, the transient topographic signal of basement exhumation has not been explored. We use model outputs to demonstrate the mobility of drainage divides in mountain ranges in response to the exhumation of basement rocks and the implications for the morphology of river catchments. The exhumation of harder rocks within a catchment reduces upstream channel steepness and erosion rates in contrast to neighboring catchments. The results are a shift in the orogen‐scale drainage divide toward the harder rocks, and the formation of range parallel longitudinal valleys as neighboring river networks capture the headwaters of catchments impacted by the harder lithology. Our model outputs provide a process explanation for the initiation of many longitudinal valleys in mountain ranges, and for the pinning of drainage divides on rocks of higher strength as seen the Central Pyrenees, Western Alps, or High Atlas.

• We investigate the impact of exhuming high strength rocks in a mountainous landscape using numerical modeling and topographic analyses • Exhumation of harder rocks causes a reduction in upstream channel steepness and erosion rate and contrast with neighboring catchments • Landscape is rebalanced by drainage divide migration toward harder lithologies and formation of longitudinal trends in river valleys

Supporting Information:
Supporting Information may be found in the online version of this article.

RESEARCH LETTER
strengths influence the erosion rate of fluvial channels by a factor of 3-5 depending on the range of lithologies (Sklar & Dietrich, 2001). Gallen (2018) shows that in the Appalachian Mountains, important drainage reorganisation can be triggered by changes in the erodibility coefficient of different rock unit types. Bernard et al. (2019) and Zondervan et al. (2020) show for the Pyrenees and High Atlas Mountains respectively, that variation of rock erodibility can counter the effect of crustal thickening or climate and force drainage divides to fix on highly resistant massifs (e.g., crystalline plutons). Studies using landscape evolution models have shown that patterns of rock uplift, climate and lithologies can lead to important changes of topography and drainage systems in active mountain ranges (Kühni & Pfiffner, 2001a;Schlunegger et al., 2001). However, understanding the landscape response to the progressive exhumation of high-strength rock, and the processes that result in divide migration and the evolution of longitudinal valleys have not been fully explored.
Here, we run a series of experiments on the impact of exhuming high-strength, basement lithologies into an active landscape, using a landscape evolution model that tracks the distribution of erosion rates. Topographic metrics, particularly the channel steepness, are extracted from the model runs to facilitate visualization of the landscape response to external changes. The time evolution of the results demonstrates the impact on drainage divides upstream of the exhumed basement rock, and on the implications for changing erosion rates and modified catchment forms.

Methods and Model Setup
In order to investigate the evolution of transient river networks in a mountainous landscape, we utilize the two-dimensional landscape evolution model Fastscape (Braun & Willett, 2013) which is available through the Xarray-Simlab package (Bovy, 2020). Erosion is simulated by river channel and hillslope processes using the stream power incision model (Howard & Kerby, 1983;Whipple & Tucker, 1999) and linear diffusion law respectively.
We coupled the Fastscape landscape model outputs with the open-source topographic analysis algorithm LSDTopotools (Mudd et al., 2014). For each of the time-steps, we analyze a series of topographic metrics that can be extracted from the Fastscape model outputs. We focus primarily on the evolution of channel steepness index ( sn k ) along river channels. The channel steepness or rate of channel slope change normalized to drainage area has been shown to reflect spatial patterns of relative rock uplift rate, precipitation or rock erodibility (Kirby & Whipple, 2012;Mudd et al., 2014;Wobus et al., 2006). We also use this metric to infer drainage divide migration supported by other topographic metrics on either side of the drainage divide such as the mean gradient and local relief (Forte & Whipple, 2018) or the hillslope relief and flow distance (Scherler & Schwanghart, 2020).
The model setup comprises an area of 450 km in length and 150 km in width ( Figure 1a). The topography initially evolves under a linear gradient of rock uplift from 0 at the lower boundary to 1 mm.yr −1 at the upper boundary over a period of 30 Myrs. The area (m) and slope (n) exponent are 0.6 and 1.5 respectively and give a concavity index ( ) of 0.4 which falls within the range of concavities (i.e., 0.35-0.6) found by studies that have analyzed the relationship between drainage area and local slope (Kirby & Whipple, 2012;Whipple & Tucker, 1999). Finally, we set a spatially uniform erodibility coefficient of 2 × 10 −6 m 0.1 .yr −1 in the range of values deduced for the sedimentary cover of mountain ranges (Gallen, 2018;Zondervan et al., 2020). When run to steady state, this setup generates an asymmetric range with the main drainage divide located in the northern part of the model, comparable to many doubly vergent mountain belts where rock uplift rates are higher on the retro-wedge due to advection of rock from the regions of maximum accretion (e.g., Taiwan and Olympic Mountains; Willett et al., 1993). This setup allows us to explore the effect of erodibility in a model where the uplift is spatially variable. As a consequence of the contrasts in rock uplift rates across the model domain, the rivers are steeper and shorter on the northern retro-wedge (Figure 1b), and are lower gradient on the southern pro-wedge (Movie S1 and S2). From 15 million years of growth to steady state, we simulate the exhumation of basement crystalline rocks by introducing a block that is 150 km long and 30 km wide with a lower erodibility coefficient (5 × 10 −7 m 0.1 .yr −1 ; Gallen, 2018). This block is located across the southern catchments with its northern boundaries at a distance of approximatively 25 km south the main drainage divide (Figure 1a).

Channel Steepness
Under steady state conditions (Figure 2a), channel steepness records the northward increase in rock uplift rates. The transient response caused by the introduction of less erodible rocks is associated with the immediate formation of two breaks in channel steepness (i.e., slope breaks G and H; Figure 2b and Figure S1c) in rivers flowing across harder rocks. From downstream to upstream, a first abrupt increase in normalized channel gradient (G; positive slope break) is observed at the southern border of the block, and a second abrupt decrease in normalized channel gradient (H; negative slope break) is located at the northern border of the block (Figure 2b and Figure 1Sc). The landscape transience is recorded by the upstream propagation of both of these slope breaks ( Figure 2b and Figures S1c-S1e). It takes approximatively 1 Myr for the upper negative slope break (H) to migrate through the whole upstream area and to reach the drainage divide (Figures S1c-S1e and Movie S2). A longer duration (∼2 Myrs) is necessary for the positive slope break (G) to migrate upstream through the harder rocks (Figures 2b and 2c, Figures S1c-S1e and Movie S2). Because of the decrease of erosion rate but continued rock uplift, river channel portions located on the harder rocks become steeper and increase in elevation as demonstrated by the channel steepness (Figures 2b and 2c) and longitudinal and transversal slope swath profile analyses ( Figures S2 and S3). Upstream river channel portions located between the harder rocks and the drainage divide also record an increase in elevation but inversely record a decrease of channel steepness (Figures 2b and 2c).
Lowering of channel steepness in the southern catchments upstream of the basement block contrasts with the surrounding unchanged channel steepness of the northern catchments (Figure 2b and Movie S2). The models indicate the vulnerability of these upper, low-gradient portions of the southern catchments to capture, particularly by the northern catchments. From 16 Myrs to about 20 Myrs (Figures 2b-2d upstream part of the southern catchment reflects a "victim catchment" (sensu Willett et al., 2014) and is progressively captured by the northern "aggressor catchment." Because rivers of the northern catchments extend southward, the slope break H is consumed at about ∼21 Myrs (Figures 2d and 2e and Movie S2). A new slope break (I) is formed when rivers of the northern catchment start to incise into the harder basement rocks (Figure 2e). From about 20 to 30 Myrs (Figures 2e and 2f) southward lengthening of the northern catchments continues until an equilibrium of channel steepness across the divide is reached, stabilizing the position of the divide.

Hillslope Erosion Rate Evolution
The distribution of erosion rate also reveals a clear response to the exhumation of harder rocks ( Figure 3). Initially, as the model runs to steady state, the distribution of erosion is determined by the regional tilt ( Figure 3a). Shortly after exhumation of the harder rocks, the upstream area of the southern catchments experiences an abrupt decrease of erosion rate, first located along the river channels but rapidly through the hillslope (Figure 3b and Figure S4). This reduction of erosion along the main river channels contrasts with the hillslope erosion, which remains relatively high due to its dependence on slope ( Figure S4); the end result is a lowering of gradients in the upstream area. In contrast, the erosion rates in the uppermost parts of the northern catchments are high, associated with southward divide migration (Figures 3c-3e). This wave of high erosion rates in the northern catchments gradually consumes the vulnerable southern catchments (Figures 3c-3e). At about 30 Myrs (Figure 3f), a new stable condition is reached and observed by a smooth gradient of erosion rate along the model and particularly along the drainage divide.

Drainage Divide Mobility
The main response generated by the important disequilibrium of both channel steepness and denudation rates (Figures 2 and 3) is the southward shift of the main drainage divide toward the vulnerable part of the southern catchment and up to the northern limit of the lower erodibility block (Figure 1a). Drainage divide migration is not instantaneous during the exhumation of harder rocks. There is a delay corresponding to BERNARD ET AL.  Figure S1 for better visualization). the time for the negative slope break (H) to propagate to the drainage divide (i.e., ∼1 Myr in this model; Figure S1 and Movie S2). The maximum drainage divide displacement rate is reached at about 18 Myrs (3 Myrs after the exhumation of the basement block at 15 Myrs) when the disequilibrium of channel steepness and erosion rate is at its greatest ( Figure 1c). After about 10 Myrs of transient landscape, the system reaches a new steady-state condition with a relatively stable drainage divide position (Figure 1c).
The prediction that drainage divides converge toward resistant crystalline massifs is supported by a number of locations. The Maladeta and Bassiès massifs in the Central Pyrenees (Bernard et al., 2019) or the Mont-Blanc massif in the Western Alps (Kühni & Pfiffner, 2001b) attract the main drainage divides.
The same pattern of disequilibrium for both the channel steepness (Movie S2) and erosion rates (Figure 4) is observed between affected and non-affected southern catchments on the left and right borders of the those containing the exhumed harder rocks. This secondary vulnerability effect results in lateral drainage capture that enables neighboring river networks to expand upstream and parallel to the block of harder rock (Figures 4b-4e). Lateral capture by neighboring rivers leads to the formation of longitudinal valleys around the margins of the harder rocks (Figure 4f). Such mechanisms may explain the occurrence of longitudinal valleys running parallel to basement massifs as observed for the upper Rhine and Rhône in the European Alps (Kühni & Pfiffner, 2001a  Myrs. The harder basement rock is initially exhumed to the surface at 15 Myrs and represented by the light pink squares with black dashed borders. Thick black lines correspond to the drainage divide between the two catchments and highlight catchment lengthening toward achieving a low contrast in erosion rates across the divide. For the first and last panels (i.e., 14.9 and 29.9 Myrs), the low contrasts in erosion rates indicate relatively stable drainage divides compared to the high erosion rate contrasts and unstable drainage divide for the intermediate panels. Initial erosion rate gradients are due to the gradient in rock uplift from south to north due to the regional tilting in the model.

Discussion
Our observations are based on topographic outputs from numerical modeling, which cannot replicate the complexity of natural settings in mountainous domains. It is clear that the model approximation used in this study misses the effects of several factors including the multiple properties of different lithologies (i.e., rock strength, fracture densities or rock weathering) Forte et al., 2016;Sklar & Dietrich, 2001;Whipple, 2002), structural faulting, climate variability or the advection of drainage divides relative to the fixed over-riding plate during active convergence (e.g., Willet et al., 2001). Fault movements can occur during the exhumation of the crystalline basement and locally perturb the landscape. Additionally, a wetter and potentially more erosive climate on the northern side of the model would have reinforced the drainage divide migration toward the basement intrusion. Finally, lateral advection of rock would further enhance the signal of southward capture as the harder block would be translated horizontally toward the central divide. However, the model simplicity where we explicitly invoke a single change in the erodibility coefficient (i.e., proxy for the rock strength), allows us to visualize the transient processes and mechanisms BERNARD ET AL.
10.1029/2020GL092210 6 of 9 happening during the exhumation of harder rocks, and hence also gain insight into the response time of these systems.
The stream power incision model is the most widely used model because of its ability to reproduce the main characteristic features of natural landscapes. However, determination of its parameters remains challenging. By testing different values of the area (m) and slope (n) exponents, we observe that the main results raised by this study (i.e., drainage divide migration and channel lateral expansion) remains consistent for an area exponent between ∼0.6 and ∼0.8 and a slope exponent between ∼1.0 and ∼1.6 ( Figures S5 and S6). For values of m and n extending beyond the ranges proposed above, the model is no longer able to propose realistic outputs. These unrealistic results are caused by the interdependency which exist between the parameters m, n, and K as previously shown by other studies (Croissant & Braun, 2014;Resentini et al., 2017;Roberts & White, 2010). The erodibility coefficient contrast between the surrounding rocks and the basement intrusion ( rocks base / K K ) plays also an important role. When the contrast ratio falls below ∼3 the basement intrusion is no longer able to capture the drainage divide and modify the channel network ( Figures S7  and S8). Similar effects occur if the basement intrusion width is too small or if the basement intrusion location is too far from the pre-existing drainage divide. When the basement width falls below ∼15 km (Figure S9) or if the basement location is further than ∼80 km from the upper boundary ( Figure S10), southern rivers are still able to cut through the basement intrusion.
The model results indicate that important localized changes in erosion rates can occur through time without varying the rate of tectonic uplift or climate. At the northern boundary of the drainage divide, the mean erosion rate is about 0.7 km.Myrs −1 during steady-state (Figure 3a). During landscape transience (Figures 3b-3e), denudation rates reach values up to 1.2 km.Myrs −1 suggesting measurable changes may solely result from divide migration. Such a change of erosion rate should be detectable with, for example, low-temperature thermochronology isotopic systems such as 4 He/ 3 He thermochronometry, which have recorded spatial variations in cooling histories for samples along fluvial valleys (Schildgen et al., 2010;Simon-Labric et al., 2014).
The increase of elevations around the stronger lithology results in a modification of the lithostatic stress distribution through the range and hence may impact the mechanical evolution of a doubly vergent wedge system (Willett et al., 1993). The result would be a change of the long-wavelength slope of the pro-wedge versus the retro-wedge. In our experiments, as the drainage divide migrates southward, the area that would correspond to the pro-wedge (southern slopes) becomes shorter and steeper and would thus encourage frontal accretion, whereas, the retro-wedge becomes longer and gentler encouraging internal thickening ( Figure 2). The change of catchment drainage area which constitutes the pro-or retro-wedge will proportionally also have an impact on the amount and rate of sediment released to the pro-and retro-foreland basins (Naylor & Sinclair, 2008). The nature of sediment eroded and resulting lithologic contents in the stratigraphic record is also expected to change as drainage divides migrate. For example, if we assume that the harder rocks in our model setup correspond to some granitoid massif, granitoid clasts will only be recorded in the pro-wedge following the timing when the basement massif is exhumed at the surface. The retro-foreland basin will record granitoid clasts after a lapsed-time corresponding to the time for the drainage divide to reach the massif (Figures 2 and 3). Drainage divide migration could therefore potentially be inferred from sediment records.

Conclusions
Based on numerical experiments, we found that the exhumation of harder rocks within mountainous drainage catchments causes a lowering of channel and hillslope gradients in the upstream area of the affected catchment. This is associated with a reduction in erosion rates, which contrasts with neighboring headwaters of other catchments that are unaffected by the exhumation of harder lithologies. This disequilibrium across catchments causes drainage divides to migrate in order to rebalance the landscape. This behavior is facilitated by an increase of denudation rate for one side of the divide, which forces the migration of the drainage divide in the direction of the harder lithologies. Additionally, neighboring divides migrate into headwaters above the exhumed harder lithology, forming longitudinal trends in river valleys. These processes are considered responsible for the pinning of the drainage divide in the central Pyrenees and for the development of longitudinal valleys such as the Rhine and Rhône rivers of the western Alps. In addition, the landscape transience related to changes of rock erodibility lasts a few millions of years as the network adapts rapidly to new conditions, and that during this time, measurable variations in erosion rates may occur without recourse to climatic or tectonic forcings.

Data Availability Statement
Data were not used, nor created for this research.