Rock strength and structural controls on fluvial erodibility: Implications for drainage divide mobility in a collisional mountain belt

18 Numerical model simulations and experiments have suggested that when migration of the main 19 drainage divide occurs in a mountain belt, it can lead to the rearrangement of river catchments, 20 rejuvenation of topography, and changes in erosion rates and sediment flux. We assess the 21 progressive mobility of the drainage divide in three lithologically and structurally distinct groups of 22 bedrock in the High Atlas (NW Africa). The geological age of bedrock and its associated tectonic 23 architecture in the mountain belt increases from east to west in the study area, allowing us to 24 track both variations in rock strength and structural configuration which influence drainage 25 mobility during erosion through an exhuming mountain belt. Collection of field derived 26


Introduction
Collisional mountains form the erosional focus of the Earth's surface.The tectonic and climatic interpretation of mountain topography and depositional stratigraphy depends on understanding the dynamics of eroding bedrock rivers.The central drainage divide of a mountain belt is the topographic boundary between river catchments draining either flank.Any movement of the drainage divide can result in the rearrangement of catchments, rejuvenation of topography, and changes in erosion rates and sediment flux (Bonnet, 2009;Giachetta et al., 2014).Bedrock erodibility is expected to play a significant role in drainage divide reorganisation since heterogeneous exhumation of weak and strong substrates can enhance and suppress erosion respectively (Giachetta et al., 2014;Gallen, 2018) and cause topographic rejuvenation, for example through the exhumation of a basement palaeosurface (Strong et al., 2019).This is especially the case in postorogenic settings where erosion is dominant over crustal thickening (Gallen, 2018;Bernard et al., 2019).However, the magnitude of erodibility variation within a mountain belt, and the mobility of the drainage divide as rivers erode through its stratigraphy are still relatively unexplored.
Collisional orogens are characterised by bedrock rivers eroding through variable rock strength and tectonic architectures, such as very strong crystalline basement, deformed strong meta-sediments, and sedimentary cover composed of weak as well as stronger strata.Recent numerical modelling studies show the complexity of incision into horizontal or gently dipping strata of varying erodibility (Forte et al., 2016;Perne et al., 2017), and also suggest steady-state denudation is more likely to develop in tilted and /or faulted strata that have been highly deformed (Forte et al., 2016).The results of these models imply drainage mobility might be prevalent during incision in sedimentary cover but not in more deformed components of a collisional mountain belt.In a field study, Gallen (2018) modelled the erosion of a hard horizontal rock layer in the Appalachians.This model predicts that a geologically instantaneous capture of the Upper Tennessee River catchment by the Lower Tennessee River occurred at 9 Ma, which has led to a shift in the drainage divide and explains observed subsequent topographic rejuvenation in the landscape visible today.In the Pyrenees, Bernard et al. (2019) show the drainage divide follows the position of high-strength, high-elevation plutons in the crystalline basement in the centre of the belt, suggesting a direct lithological control on the position of the central drainage divide.These field studies demonstrate the drainage divide can be mobile in near-horizontal sedimentary stratigraphy of mountain belts and is likely to move to the centre of highly resistant plutons as they get exhumed in the axis of a collisional mountain belt.
However, both studies also demonstrate the challenge of documenting drainage mobility, which is often only recognisable in instantaneous capture events which leave pervasive topographic evidence.In addition, whilst numerical simulation studies are instrumental in predicting long-term processes in geomorphology which are hard to derive from observations of modern day landscapes alone, these involve simplifications and it is challenging to model the full complexity of lithologies and structural geology of a collisional mountain belt.Furthermore, the fluvial erodibility of rock, which depends on rock mass strength as well as jointing and weathering, is challenging to measure quantitatively (Bursztyn et al., 2015) and consequently the range of erodibility inputs in modelling studies vary widely (Roy et al., 2015;Forte et al., 2016;Yanites et al., 2017).
Consequently, while theory and numerical models predict a role for bedrock erodibility in driving drainage divide mobility, where and when this is expected to occur in the evolution of mountain belts needs to be understood from a range of settings with varying rock strength conditions.The challenges of measuring and characterising bedrock erodibility (Bursztyn et al., 2015) and drainage divide dynamics in field settings (Willett et al., 2014;Forte and Whipple, 2018) remain key problems.
Addressing this challenge requires a dataset on rock erodibility that is combined with topographic measures of drainage divide mobility in a mountain belt where erosion is dominant over tectonic advection, and where rivers erode through a lithologically variable landscape.In this study, we aim to derive data and present it in a form that is useful for both field geologists and numerical modellers alike.We focus our study on the central High Atlas Mountains (Fig. 1), where different stages of orogenic landscape evolution, from erosion through gently dipping sedimentary cover to exhumed crystalline basement occur along the length of the main drainage divide, forming a natural experiment.Furthermore, the continental inland setting, low weathering rates associated with an arid climate and low rates of tectonic activity make this post-orogenic mountain range an ideal setting to study erodibility-induced drainage divide mobility.In this study, we quantify the magnitude of variation in erodibility between rock types by: (i) extracting the normalised river channel steepness of rock units from a digital elevation model; (ii) collecting mechanical measurements of rock strength in the field, and (iii) quantifying the mobility of the drainage divide using topographic metrics.Our measures of rock erodibility and stratigraphic orientation are then used to quantitatively assess when and where a drainage divide will move in response to erodibility 104 Headwater channel locations used for calculating drainage divide mobility start at a reference drainage area of 1 km2, marked as white dots (see methods 3.3).The drainage divide is marked on as a black line and is segmented into five equal lengths (Fig. 5).b) Cross-sections based on seismic sections modified from Errarhaoui (1998); Teixell et al. (2003).These cross-sections are vertically exaggerated by a factor of 1.5.The location of the drainage divide is marked on as black markers on the structural cross sections.
variation in the evolution of a collisional mountain belt in which rivers erode first through sedimentary cover, secondly through meta-sediments and finally through crystalline basement.
Results from this study have implications for understanding long-term trends in the geomorphology and erosion of mountain belts owing to the rock types and their structural configuration, and quantitatively constrain the magnitude of contrasts in erodibility within a mountain belt which enables integration of observation and modelling studies.

Study area
The study area is located in the centre of the High Atlas Mountains, focussed on the 250-300 km length of the drainage divide bounded by the Ouarzazate Basin to the south and the Haouz Basin in the north (Fig. 1).The drainage divide strikes along the structural grain of the mountain belt, inherited by the configuration of a pre-existing rift (Babault et al., 2012).The age of exposed bedrock and its associated tectonic architecture in the mountain belt increases from east to west in the study area (Fig. 1).The along strike trends allow us to use location as a proxy for time in understanding the mobility of the main drainage divide as river erosion exposes sedimentary cover to metamorphic sedimentary rock, to crystalline basement.The sedimentary cover, older meta-sediments and underlying crystalline basement have distinct lithologies with varying hardness and structural weaknesses owing to the tectonic architecture of the mountain belt (Table 1).Thus, each chronostratigraphic package has its own potential for drainage divide mobility which we investigate in this study.
From east to west the drainage divide crosses first a landscape dominated by gently dipping beds of Mesozoic carbonates and clastic sedimentary rock (Table 1) configured into low-amplitude longwavelength folds punctuated by spaced-out thrusts (Fig 1,Fig 2a) at an elevation of 2600 to 3000 m.
In this length of the divide, the Dades catchment lies to the south and incises into a long-wavelength syncline exhibiting gently dipping strata, whereas directly north of the divide is a concentration of folding and thrust offset strata dipping at higher angles to the surface (Fig. 1b).Next the drainage divide decreases in elevation to ~2500 m, as it crosses a zone of dipping strata of faulted metasediments (Table 1, Fig. 2b) in the centre of the study area (Fig. 1).Finally, in the west, the divide rises back up to elevations > 2800 m as it runs over the edge of the exhumed crystalline basement consisting of igneous rocks with vertical or sub-vertical fracturing (Fig. 1, Table 1, Fig 2c-d).
The inland setting of the central High Atlas (800-1000 km from the coast) results in an absence of eustatic base level control on drainage development.Though localised base level fluctuations are likely to have occurred, studies of basin fill sediments and river profiles suggest that such fluctuations were small to negligible over the Plio-Quaternary (Boulton et al., 2014;Boulton et al., 2019).The High Atlas is in a post-orogenic state, with long-term isostatic rock uplift rates of 0.17-0.22mm yr -1 since 15 Ma related to lithospheric thinning (Babault et al., 2008).The lack of recent tectonic deformation is apparent in the undeformed continuous Quaternary river terraces forming parallel river long profiles throughout the fold-thrust belt and thrust front (Stokes et al., 2017).
Therefore, unlike in active mountain belts where divide migration is expected in response to changes in the uplift field (Willett et al., 2001), the Plio-Quaternary divide mobility of the High Atlas is not expected to reflect tectonic advection.Recorded glacial features related to Pleistocene glaciation establish that the snowline was at c. 3300 m in the High Atlas (Hughes et al., 2004).Since the main drainage divide varies in elevation between 2500-3000 m it has not been affected significantly by glacial activity.Thus, for our purposes the main control on the Plio-Quaternary evolution of the High Atlas river network is the incision through the inherited tectonic architecture of lithological units and their contrasting strengths.Drainage development in the High Atlas is considered to be primarily the product of the exhumation of structurally distributed lithologies with different hardness, controlling where river terraces develop (Stokes et al., 2017) and affecting the occurrence of diffusive and advective slopes (Mather and Stokes, 2018).The High Atlas is set in a semi-arid climate, which means that the effect of weathering on rock erodibility is expected to be low.Thus, intact rock strength measurements are likely to reflect the effective bedrock strength.

Methods
We measure two proxies of fluvial bedrock erodibility to constrain the contrast between lithological units.These data allow us to compare our field observations of drainage divide mobility with numerical models of river erosion through variable lithology.We systematically: 1) derive the normalised river channel steepness index from a digital elevation model (DEM) as a measure of a river channel's power to erode rock; 2) record compressive rock strength using a Schmidt hammer in the field as another measure of fluvial erodibility, and 3) extract topographic metrics of drainage divide mobility from the DEM along the length of the main drainage divide.

River profile analysis, rock type and fluvial erodibility
Rock type influences the river network by affecting the ability of rivers to erode into bedrock, determined by the fluvial erodibility.Numerical models that describe river erosion through bedrock often use the stream power model, in which the erosion rate at any particular point in a bedrock river channel is defined by: where K is an erodibility constant which depends on the rock-type over which the river channel flows as well as the climatic setting, A is upstream drainage area and S is local channel gradient, and m and n are constants that depend on basin hydrology, channel geometry, and erosion processes (Whipple and Tucker, 1999).Any change in the erodibility of rock type exposed will force the river to adjust its stream power by changing river channel slopes and thus can cause divide mobility as river networks respond.The use of this formula in field studies proves problematic because of the lack of data and robust methodology to determine the erodibility constant, K.However, equation 1 can be written as which may be recognised as a form of the empirical power law scaling local channel gradient (S) and drainage area (A): =    − (Eq. 3) where S is the local channel slope (dimensionless), ksn is normalised steepness index (m 2θ ) and  is the concavity index (dimensionless).We assume that erosion is proportional to specific stream power and inversely proportional to bedrock erodibility, so that n > 1 (Perne et al., 2017), and ksn is inversely proportional to erodibility, K: Unlike K, the ksn of river channels can be readily determined from digital elevation models (DEMs) (e.g.Boulton et al., 2014;Gallen, 2018;Bernard et al., 2019).We therefore use ksn as a measure of fluvial erodibility of geological units over which river channels flow in the High Atlas, provided that spatial variability in long term rock uplift is low, which is likely, for reasons outlined above (Section 2).Consequently, a factor of difference in ksn between geological units will be an estimate of the factor of difference in their fluvial erodibility, K.We normalise ksn values to the most erodible rocktype (see section 4.1), and then use Eq. 4 to convert these normalised average ksn values to normalised K values.To do this, we calculate normalised K for three feasible values of n: n=1, n=2 and n=4.Perne et al. (2017) show that to obtain a stream-power-model generated river profile in which the slopes of river channels are steeper in low erodibility, strong rock, n must be > 1. Lague (2014) demonstrate that, fully calibrated with slope information from field locations, n ~ 2 in most cases, and this could be as high as n ~ 4. Calculating normalised K values from the ksn data enables a quantification of contrasts in erodibility K between rock-types in the High Atlas and their control on changes in river erosion rates and consequent mobility of the drainage divide.
We performed ksn analysis using the Topographic Analysis Kit (Forte and Whipple, 2019), a series of MATLAB functions based on TopoToolbox (Schwanghart and Scherler, 2014) which uses the chi approach to calculate ksn smoothed over 1000 m segments (Perron and Royden, 2013).Tests using the criteria proposed by Perron and Royden (2013) show that  is ~ 0.45 for a range of catchments covering the extent of the drainage divide in the central High Atlas, which is the same value as used for the High Atlas by Boulton et al. (2014).This gives ksn units of m 0.9 .We use the ALOS Digital Elevation Model with a resolution of 30 m from the Japan Aerospace Exploration Agency (http://www.eorc.jaxa.jp/ALOS/en/aw3d30/index.htm) to perform the analysis following recommendations of Boulton and Stokes (2018).Finally, we extract the average and standard deviation ksn of each rock type using the digitised 1: 1 000 000 geological map of Morocco (Saadi et al., 1985, see Supplementary Fig. S1).Formations are grouped into lithological units that collectively comprise three groups: Mesozoic sediments, Palaeozoic meta-sediments and crystalline basement over which the drainage divide crosses from east to west (Fig. 1).For each of those packages we group lithological units (Table 1): (i) gently folded massive marine platform limestone; (ii) interbedded red beds consisting of mud and siltstone; (iii) interbedded sand-siltstones and conglomerates; (iv) folded and faulted meta-sedimentary rocks such as schist and shale, and (v) faulted massive outcrops of igneous units.The contrasts in ksn of river channels flowing over rocktypes (Supplementary Fig. S1) includes the effects of both rock strength and structural weaknesses on the resulting erodibility of bedrock.

Mechanical rock strength measurements, rock type and erodibility
A further approach to determine the fluvial erodibility of rock-types is measuring their mechanical strength.The classic experiment by Sklar and Dietrich (2001) showed that erodibility of rock, K, in the stream power model for bedrock river erosion is related to the inverse square of tensile strength (  2 ).Thus, a measure of tensile strength and the difference between units enables the calculation of contrasts in erodibility, K.However, measurements of tensile strength cannot be achieved easily in the field, and the most commonly applied method of rock strength assessment in geomorphology is the Schmidt hammer because of its portability (Goudie, 2006).The Schmidt hammer records the rebound distance of a weighted spring that strikes the rock surface and uses a range from 10 -100.
With the Schmidt hammer, measurements of in situ uniaxial compressive rock strength in the landscape can be taken in large quantities, making it a very versatile instrument in landscape studies.
The higher the rebound value, the higher the elastic strength of the rock, which is a measure of the uniaxial compressive strength (UCS) of a rock.Tensile strength (TS) represents the resistance to sediment impacts on the riverbed and its use has been advocated based on the stronger correlation of fluvial metrics with tensile strength than compressive strength (Bursztyn et al., 2015).But since compressive and tensile strength are correlated (UCS ~ 10*TS; Kahraman et al., 2012;Nazir et al., 2013), K is effectively also proportional to the inverse square of UCS:  ∝ 1  2 (Eq.5) Thus, similar to the approach in 3.1, we normalise UCS values to the most erodible rock-type, and then use Eq. 5 to convert these normalised average UCS values to normalised K values.To produce a map of rock strength we digitised the 1:1 000 000 geological map of Morocco (Saadi et al., 1985) for the study area and assigned an average rock rebound value to each stratigraphic unit.Typically, ten to twenty Schmidt hammer measurements were taken at each location, totalling 690 readings throughout the central High Atlas with up to 132 readings per geological unit (Supplementary materials).Where rock strength measurements for geological units from field measurements were absent, we supplemented this with an existing database of Schmidt rebound values for lithological units reported from around the world (Goudie, 2006).This database contains Schmidt rebound measurements for 111 lithological units, in mostly arid environments similar to the High Atlas.For each stratigraphic unit without field measurements the average of values reported in the literature for the lithology of that unit is used.Units are combined in the same lithological groups as the ksn data.The standard deviation of Schmidt hammer measurements and the range of values reported in the literature reflect the variation of rock strengths within the lithological units (Supplementary material Table 1, Figure S2).Mean and standard deviation of Schmidt hammer rebound values (SHV) are then converted to UCS using the conversion which was derived by Katz et al. (2000) ( = 2.21 (0.07) ) for a range of carbonate rocks, sandstone, marble and igneous rocks with UCS values similar to those found in the study area.In contrast to the ksn approach, the UCS of lithological units only takes into account the internal rock mass strength of bedrock, and thus does not account for structural weaknesses imposed by discontinuities such as bedding or jointing.Consequently, estimating variation in erodibility between rock-types using this approach may be expected to underestimate the erodibility of folded and faulted shales and schists and jointed igneous rock (Table 1).

Topographic analysis of the drainage divide
To determine the mobility of the drainage divide in response to erodibility variation and its structural configuration, we perform topographic analysis of the main drainage divide along the  Forte and Whipple (2018).Where channels are steeper on one side of the divide compared to the other, divide motion will be in the direction of low gradients where erosion is lower.B) Similarly, erosion will progress towards those channels with higher headwater channel elevations.
study length.The mobility and potential direction of movement of the drainage divide depends on erosion rates either side, and whilst the chi method of mapping drainage divide instabilities (Willett et al., 2014) is used widely, Forte and Whipple (2018) demonstrated that this method is especially problematic when integrating across multiple lithologies with different strengths.Alternatively, Forte and Whipple (2018) coined the term Gilbert Metrics, based on the assumption that contrasting stream gradients either side of a divide will drive divide migration.Where divides are bounded by different channel gradients either side they are mobile, with higher erosion rates on the steeper sides leading to migration towards the side with lower channel gradients (Fig. 3).The topographic proxies for erosion rates across divides defined by Forte and Whipple (2018) are differences in headwater channel elevation, local headwater hillslope gradient and local headwater relief, which have proven useful for interpreting the relative mobility of catchment divides (Forte and Whipple, 2018).Headwater channel elevation, local headwater hillslope gradient and local headwater relief were extracted from the ALOS DEM using methods outlined in Forte and Whipple (2018), at a reference drainage area 1 km 2 .While the upstream area required to initiate channel flow varies on settings (Clubb et al., 2014), this drainage area is the critical threshold area downstream of which stream-flow-dominated fluvial channels are expected to dominate over debris-flow-dominated colluvial channels in most settings (Wobus et al., 2006).In a numerical landscape evolution model, Forte and Whipple (2018) test the applicability of the Gilbert metrics in a setting with erodibility contrasts across a drainage divide, and conclude that divide migration rate approximates a linear relationship with cross-divide differences in erosion rates and all three Gilbert metrics.To assess the mobility of the drainage divide we define five equal lengths of about 50 km each where boundaries align with those of the three lithological packages described earlier (Fig. 1), and in each the mean, standard deviation, standard error and bootstrap confidence interval of headwater channel values are calculated either side of the drainage divide.2, modified from Goudie (2006).Erodibility classification is based on rock strength (Table 2).

Ksn and fluvial erodibility
From east to west, the drainage divide runs over moderately resistant to erodible and very erodible sedimentary cover, moderately resistant meta-sediments, and resistant crystalline basement (Fig. 1).
Red beds, other clastic sediments and limestones in the Mesozoic cover have mean ksn values of 37, 59 and 74 m 0.9 with standard deviations of 25, 42 and 49 m 0.9 respectively (Fig. 4c).Palaeozoic metasediments average at 85 m 0.9 with a standard deviation of 41 m 0.9 , whilst the crystalline basement has the highest ksn at 146 with a standard deviation of 52 m 0.9 (Fig. 4c).Fluvial erodibility K values based on Eq. 4 and ksn, normalised to the Mesozoic red beds (the weakest lithology: section 3.2), vary by a factor of four to fifteen (Fig. 4d).The normalised K value for red beds is 1.0, with other clastics 0.40 or 0.63 and limestone in the Mesozoic cover at 0.25 or 0.50, depending whether n = 2 or n = 1 (Fig. 4d).Normalised K for Palaeozoic meta-sediments is 0.19 or 0.44 whilst crystalline basement has the lowest normalised K at 0.066 or 0.26, about four to fifteen times less erodible than the Mesozoic red beds.

UCS and fluvial erodibility
Similar to the ksn of rock units, the fluvial erodibility of bedrock along the length of the drainage divide based on UCS varies from alternating very erodible to moderately resistant rock in the sedimentary cover, to moderately resistant in the meta-sediments and resistant in the crystalline basement.The Mesozoic sedimentary cover has the largest range in UCS of the three chronostratigraphic packages, with red beds, other clastic sediments and limestones averaging 12, 36 and 57 MPa respectively, with standard deviations of 7, 12 and 18 MPa (Fig. 4a).UCS values for the Palaeozoic meta-sediment averages at 55 MPa with a standard deviation of 19 MPa, similar to the Mesozoic limestones (Fig. 4a).The igneous rocks of the crystalline basement have the highest UCS values with a mean and standard deviation of 101 and 20 MPa.Grouping of rocks into five classes is based on a table from Goudie (2006) which is based on Schmidt hammer values and other measures of rock strength (Table 2), and classes are marked on the UCS graph (Fig. 4a).Fluvial erodibility K values based on Eq. 5 and UCS, normalised to the weakest unit, vary by two orders of magnitude (Fig. 4b).The normalised K value for red beds is 1.0, with other clastic and limestone in the Mesozoic cover at 0.1 and 0.05 respectively.Normalised K for Palaeozoic meta-sediments is 0.05 whilst crystalline basement has the lowest normalised K at 0.015, about two orders of magnitude less erodible than the Mesozoic red beds.Categories of erodibility suggested based on the lithological grouping into rock strengths (Table 2) and their respective normalised K values based on UCS are defined in Fig. 4b.

Drainage divide mobility
Normalised cross-divide delta values of the Gilbert metrics following Forte and Whipple (2018) show the magnitude and direction of mobility varies along the divide and that the divide is stable in the central length of the divide (Fig. 5).Elevation of the drainage divide is high in both the east and west, around 2600-3000 m in DIV 1 and DIV 4-5, and lowest in the middle length at ~ 2500 m (Fig. 5a).This low elevation length (DIV 2-3) of the divide in the centre of the study area has equal headwater hillslope gradients and relief on both sides of the divide (Fig. 5b-c) suggesting a stable drainage divide in the meta-sedimentary rocks (Fig. 1, 5d).From west to east, headwater hillslope gradients and relief decrease from 0.6 (35°) and 500 m in the crystalline basement of DIV 1 to 0.3 (15°) and 220 m in the Mesozoic sedimentary cover of DIV 5. Values of headwater channel elevation, hillslope gradient and relief differ either side of the drainage divide in the eastern (DIV 1) and western (DIV 4-5) lengths (Fig 5b-c).In DIV 1, the Gilbert metrics indicate northward movement of the divide (Fig. 5d) towards the centre of the exposed crystalline basement.This mobility is based on headwater values which are on average 169 m higher in elevation to the north, with lower hillslope gradients and relief, which differ by 0.06 (3.4°) and 50 m of relief (10%) across the divide (Fig. 5a-c).In DIV 4 and 5 in the Mesozoic cover the divide is significantly mobile towards the south (Fig. 5d), with headwater channel elevation contrasts of 96-113 m, a difference in hillslope gradients of 0.07-0.09(4-5°) and relief 44-55 m (14-30%) (Fig. 5a-c).

Rock strength, ksn and fluvial erodibility
The theory and empirical relationships outlined in sections 3.1 -3.2 predict that fluvial erodibility, K relates to UCS in an inverse square (Eq.5) and to ksn in an inverse relationship in which the power depends on the value of n (Eq.4).We find a linear relationship between UCS and ksn, which suggests that n = 2, such that ksn and UCS scale according to a linear relationship.However, visual inspection of figures 4c and 4d suggests that normalised K values calculated from ksn data where we assume n = 4 are more similar to the normalised K values derived from UCS data.However, since Lague (2014) found that n ~ 2 is most commonly observed in the field, and our regression of UCS and ksn data yields a strong linear fit consistent with n = 2, we suggest that n ~ 2 and that any difference between UCS and ksn-derived normalised K values is due to other effects.For example, UCS does not explicitly include other factors influencing bedrock erodibility including the degree of weathering and structural discontinuities (Table 1), which especially through zones of deformation, will lead to more rapid erosion of even hard rocks (high UCS) by fluvial systems.For example, Anton et al. (2015) and Baynes et al. (2015) showed that canyons can be created by extreme flood events in basalt and granite respectively, where the presence of discontinuities enables rapid erosion through fluvial plucking and block topple.On the other hand, whilst the stream power model of bedrock river erosion only accounts for changes in river channel slope, field studies show that rock strength correlates with channel width (Allen et al., 2013), as well as valley width (Schanz and Montgomery, 2016) and can influence the efficiency of river bed load in eroding underlying bedrock (Brocard and van der Beek, 2006).Furthermore, there can be a dampening of ksn value variations across lithological boundaries as sections of river with weak bedrock downstream of river reaches with hard bedrock can be armoured with blocks (e.g.Thaler and Covington, 2016).Based on the lithological effects on river channel and valley morphology demonstrated by these field studies, using ksn as a measure of rock erodibility in the stream power model of bedrock river incision likely underestimates the effect of lithology on river erosion.
Thus, calculated through UCS measurements, K is expected to vary by two orders of magnitude (Fig. 4c), whereas if using ksn, K is expected to vary by one order of magnitude only (Fig. 4d).The lack of constraints of K in natural settings have led to numerical modelling studies varying widely in the range of erodibilities, using from one (Forte et al., 2016), two (Yanites et al., 2017), to three orders of magnitude difference between rock-types (Roy et al., 2015).Such a range of inputs is often based on the experimental relationship between intact rock strength and erosion in the classic abrasion mill experiment done by Sklar and Dietrich (2001), and the five orders of magnitude difference in K derived through the early work of Stock and Montgomery (1999).The latter forward-modelled river paleo-profiles, constrained by bedrock strath terraces and basaltic layers, to presently observed profiles for a range of locations worldwide.Their values of K range from 10 -2 to 10 -4 m 0.2 yr -1 in the mudstones of humid continental Japan, to 10 -6 -10 -7 m 0.2 yr -1 in the subtropical granite landscape of Australia (Stock and Montgomery, 1999).An issue with the experimental approach by Sklar and Dietrich (2001) is that it does not include the effects of weathering and jointing of rock in natural landscapes and their influence on fluvial erodibility.In the results from Stock and Montgomery (1999) the influence of rock-type is difficult to isolate from the variation in climatic setting owing to the spread of study locations.Thus, whilst these early studies give some first order estimates of possible absolute values of fluvial erodibility and the relationship between rock strength and erosion rates, our results constrain more fully the contrasts in fluvial erodibility between rock types which may be expected within a mountain belt.
Next to rock strength control on ksn, orographic enhancement of precipitation with greater precipitation on the northern and western sides of the drainage divide may lead to in a decrease of ksn values from south to north or east to west (Supplementary Fig. S3).However, there is no evidence to suggest a significant difference in ksn values between north and southern or east and western portions of the High Atlas (Supplementary Fig. S4).

Drainage reorganisation in sedimentary cover
The Gilbert metrics indicate drainage divide mobility where it crosses the sedimentary cover in the east (Fig. 5d).Here, the Mesozoic sedimentary cover is gently deformed, resulting in gently dipping strata punctuated by widely spaced thrusts and folds (Fig 1b,Fig 2a).For example, the Dades river catchment south of the main drainage divide (Fig. 6a) incises into a long-wavelength syncline of The erosion and exposure of a hard limestone layer underneath a soft red bed layer in the Mesozoic low amplitude syncline leads to a divide mobility towards the southern catchment (Fig. 5).Panel a) shows the landscape at time t0 with the location of the divide between southern catchments with underlying sub-horizontal strata and northerly catchments with steeply dipping strata.b,c) As the rivers erode through the sub-horizontal strata of the southern catchment the ksn of the surface exposed needs to adjust to the new lithological strength, steepening transiently as it incises into harder limestone bedrock.During the response, erosion rates in the catchment drop before re-establishing equilibrium with uplift rates.On the other hand the northerly catchments do not experience such a change in ksn and experience steady rates.d,e) Graphs show the change in catchment-wide erosion rates either side of the drainage divide from t0 to t1 when the lithological contact gets exposed in the southern catchment to time t2 when erosion rates have reached steady state.The graph for the southern catchments shows the results from landscape evolution modelling of the exposure of a contact between sub-horizontal soft stratigraphy on top of a hard layer (Forte et al., 2016).The resultant change in erosion rates across the divide explains the migration of the drainage divide.The response depends on the factor of difference in erodibility, K, between units.Panel f) shows how the erosion of the lithological contact in the southern catchments leads to drainage divide migration to the south over the time of the transient response.
slightly dipping strata composed of weak continental red beds and hard limestones, whereas directly north of the divide folding and thrusting is more closely spaced resulting in strata dipping at higher angles to the surface (Fig. 6a).Our results show that ksn is correlated to rock strength (Fig. 4e), and in the hard limestone ksn is higher than in the weak red beds (Fig. 4a,c).Whereas erosion through nearvertical strata north of the divide result in near-stable ksn values through time (Fig. 6b), the horizontal stratigraphy of the Dades river catchment to the south of the divide means ksn values need to change through time to return to stable erosion rates (Fig. 6c) which equal rock uplift rates.
Consequently, there is a period of transience when ksn values adjust to the change in bedrock erodibility that occurred when erosion of soft red beds exposed hard limestone along the majority of the river catchment (Fig. 1a, 6f).This transition period explains the southwards movement of the divide (Fig. 5d), as erosion rates stay more or less stable in the north (Fig. 6d) whilst transient ksn values cause a temporary decrease in erosion rates in the southern river catchments (Fig. 6e).The stratigraphic effect on river erosion rates presented in Fig. 6e was first demonstrated in a numerical modelling study of river erosion through layered stratigraphy by Forte et al. (2016), who show that for strata of variable erodibility dipping 5° or less the overall erosion rate of the landscape is expected to change by several factors during incision.Here, we show that such an effect can lead to migration of the drainage divide in the gently deformed sedimentary cover in collisional mountain settings (Fig. 6f), where exhumation of hard and soft strata is isolated or offset by faults.
Simulations by Forte et al. (2016) show the exhumation of a stratigraphic contact with a factor of two to ten difference in erodibility (Fig. 6e) could take 2 -9 Ma to re-equilibrate in an area of 800 km 2 with a relatively high amount of precipitation (1 m yr -1 ).Therefore, drainage divide migration in response to the incision through soft red beds to hard limestones in the upper half of the 1500 km 2 Dades catchment, representing a change in erodibility by a factor of 2 -20 (Fig. 4), where rainfall is on the order of 0.1 -0.5 m yr -1 , is expected to persist on a timescale of 10 6 -10 7 yrs.

Divide migration driven by crystalline basement exhumation
Where crystalline basement is exposed we find the position of the drainage divide is shifting towards the centre of this exposed resistant rock in the High Atlas (Fig. 5).In a numerical simulation, Giachetta et al. (2014) found that when they imposed an erodibility gradient across a drainage divide, which is representative of the exhumation of a crystalline basement such as in the High Atlas, the drainage divide responded by moving towards the side of lower erodibility over a timescale on the order of 10 6 yrs.Bonnet (2009) proposed that such a shift of the drainage divide is accompanied by a split of catchments, creating more and smaller catchments on the side of lower erodibility.Giachetta et al. (2014) also show that on the other side, a growth of larger catchments occurs.
Similarly, Bernard et al. (2019) found that the drainage divide in the Pyrenees follows the position of strong plutons.This implies that our results show a transient stage of drainage divide migration in response to exhumation of crystalline basement, where today's drainage divide at the edge of the crystalline basement is expected to be stable in the centre of the strong crystalline basement, or might even continue reorganising within the basement to follow the exhumation of resistant plutons.Giachetta et al. (2014) used two orders of magnitude difference in erodibility values to model this effect, and here we show divide mobility can be driven by exhumation of basement that is only a factor of two less erodible than the overlying meta-sedimentary rock if calculated through ksn (Fig. 4d), and a factor of three less erodible if calculated through UCS (Fig. 4b).

Lithologically induced drainage divide mobility during the long term erosion of a collisional mountain belt
The combination of estimations of contrasts in fluvial erodibility of rock types (Fig. 4), geomorphic indicators of drainage mobility (Fig. 5), and considerations of their structural configuration in the rift such as in the High Atlas (Babault et al., 2012).We purport that changes in erosion rates as rivers incise through strata of different erodibility will drive drainage reorganisation in collisional mountain belts, where layers are close to horizontal and only gently deformed (Fig. 6, 7).This is because where strata are deformed gently and offset by faults, local exhumation of contacts between soft and hard rock leads to changes in erosion rates between catchments (Fig. 6).Consequent changes in erosion rates across the drainage divide will lead to the migration of the drainage divide (Fig. 7a) as illustrated by the mobility of the drainage divide in the High Atlas (Fig 5,6), which could lead to steady divide migration or instantaneous capture of catchments, such as shown in the Appalachians (Gallen, 2018).The effect of rock type on drainage reorganisation will be strongest in early phases of collisional mountain building, before the sedimentary cover erodes in the centre of the belt.In later stages, minor reorganisation and capture could still occur closer to the thrust front on the margins of the mountain belt where Mesozoic sedimentary strata are present.When deformed meta-sediments become exhumed, the increase in dip and deformation of strata leads to more stable erosion as rivers incise (Forte et al., 2016, Fig. 6), resulting in stable drainage divides as witnessed in the middle of the study area (Fig 5,7b).When crystalline basement gets exhumed, the drainage divide will migrate into the centre of highly resistant rocks resulting in more drainage divide migration (Fig 5,7c).The migration of the main drainage divide in a mountain belt has been shown to lead to reorganisation of river catchments (Bonnet, 2009;Giachetta et al., 2014), imposing new boundary conditions on river channels which change gradients and sediment loads (Forte et al., 2015) and the ensuing response can result in a cascading effect, impacting geomorphic and stratigraphic systems for millions to tens of millions of years (Beeson et al., 2017).The effects of lithologically-induced drainage migration here described for a post-orogenic belt could be more complex in an active mountain belt setting.

Conclusions
This study shows that in a collisional mountain belt, the drainage divide will be mobile in response to changes in erosion rates of rivers incising into gently dipping and deformed strata of contrasting erodibility in the sedimentary cover, and in response to the exhumation of strong crystalline basement.A combination of rock mass strength measurements and ksn derived from a digital elevation model constrain the contrasts in fluvial erodibility exhibited in the High Atlas to between a factor of 4 calculated through ksn and two orders of magnitude calculated through UCS.In the stage of a collisional mountain belt during which rivers incise through highly deformed meta-sedimentary units of intermediate erodibility, rock-type induced boundary conditions affect river networks the least.Based on our values of erodibility contrast and previous numerical models we estimate the timescale of adjustment in response to changes in erodibility of exposed bedrock to be on the order of 10 6 -10 7 yrs.Our results demonstrate that the mobility of the drainage divide in a collisional mountain belt can be driven by rock erodibility variation alone, which has implications for the perception of autogenic dynamism of drainage networks and fluvial erosion in collisional mountain belts, and the interpretation of their geomorphology and downstream stratigraphy.

Figure 1
Figure 1 -a) Geological map of the central High Atlas, showing the distribution of Mesozoic sedimentary cover consisting of limestones (blue) and continental clastic sedimentary rocks (pink), Palaeozoic meta-sediments (orange) and crystalline basement (red).The chain is flanked by sedimentary basins: the Ouarzazate Basin to the south and the Haouz Basin to the North.Three cross-sections labelled A to C show the locations of cross-sections in b).Headwater channel locations used for calculating drainage divide mobility start at a reference drainage area of 1 km2, marked as white dots (see methods 3.3).The drainage divide is marked on as a black line and is segmented into five equal lengths (Fig.5).b) Cross-sections based on seismic sections modified fromErrarhaoui (1998);Teixell et al. (2003).These cross-sections are vertically exaggerated by a factor of 1.5.The location of the drainage divide is marked on as black markers on the structural cross sections.

Figure 2 -
Figure 2 -field photos showing the lithologies and stratigraphic configuration of units within the three lithostratigraphic packages of the Atlas.a) In the Mesozoic weathering-resistant competent limestone beds dip at subhorizontal angles.Within the interlayered limestones and weaker clastic rocks, variations in rock strength and weathering exist.b) In the Palaeozoic meta-sediments steeply dipping, sub-vertical strata are common.Road on the flood plain for scale.c) Massive crystalline basement exposed in is sometimes found jointed in vertical plains or d) as dome-shaped or massive igneous outcrops with vertical or sheet fracturing.House for scale at the bottom.

Figure 3 -
Figure3-principles of the Gilbert metrics for drainage divide stability, modified afterForte and Whipple (2018).Where channels are steeper on one side of the divide compared to the other, divide motion will be in the direction of low gradients where erosion is lower.B) Similarly, erosion will progress towards those channels with higher headwater channel elevations.

Figure 4 -
Figure 4 -Average and standard deviation of a) uniaxial compressive strength b) ksn for every chronolithological unit.The number of values analysed for each lithology is displayed above the plots.Colours in graphs are in accordance with Fig. 1.Sources of compressive strength data: 1)Gokceoglu and Aksoy (2000) 2)Goudie (2006) and references cited therein 3)Goktan and Gunes (2005) 4)Karakus et al. (2005) 5)Pye et al. (1986) 6)Kahraman et al. (2002).Fluvial erodibility K for each chronolithological unit normalised against the weakest rock type derived from relative UCS values (c) and relative ksn values (d).In d), values of normalised K derived from ksn depend on the value of variable n in the stream power equation (Eq.4) e) Linear regression of values of ksn and UCS for every geological unit (geological map 1:100 000) with a 95 % confidence envelope in dashed grey lines and points coloured by chronolithological membership.R is the correlation coefficient.Classification of rock strength is based on Schmidt hammer values and other measures of rock strength summarised in Table2, modified fromGoudie (2006).Erodibility classification is based on rock strength (Table2).

Figure 5 -
Figure 5 -Gilbert metrics for each length of the drainage divide (Fig. 1) from west to east.a-c) headwater channel mean and standard error elevation, gradient and relief respectively.d) Potential direction of divide migration indicated by the direction of mobility of the divide from the Gilbert metrics.Values are standardized to show the direction of mobility.Bars show the standard deviation and shaded boxes show bootstrap confidence intervals.

Figure 6 :
Figure6: The erosion and exposure of a hard limestone layer underneath a soft red bed layer in the Mesozoic low amplitude syncline leads to a divide mobility towards the southern catchment (Fig.5).Panel a) shows the landscape at time t0 with the location of the divide between southern catchments with underlying sub-horizontal strata and northerly catchments with steeply dipping strata.b,c) As the rivers erode through the sub-horizontal strata of the southern catchment the ksn of the surface exposed needs to adjust to the new lithological strength, steepening transiently as it incises into harder limestone bedrock.During the response, erosion rates in the catchment drop before re-establishing equilibrium with uplift rates.On the other hand the northerly catchments do not experience such a change in ksn and experience steady rates.d,e) Graphs show the change in catchment-wide erosion rates either side of the drainage divide from t0 to t1 when the lithological contact gets exposed in the southern catchment to time t2 when erosion rates have reached steady state.The graph for the southern catchments shows the results from landscape evolution modelling of the exposure of a contact between sub-horizontal soft stratigraphy on top of a hard layer(Forte et al., 2016).The resultant change in erosion rates across the divide explains the migration of the drainage divide.The response depends on the factor of difference in erodibility, K, between units.Panel f) shows how the erosion of the lithological contact in the southern catchments leads to drainage divide migration to the south over the time of the transient response.

Figure 7 -
Figure 7 -Conceptual model of the development of a collisional mountain belt and the behaviour of the central drainage divide in response to exhumation of lithostratigraphic units

Table 1 -
chronostratigraphic packages labelled on Fig 1 with details of lithologies and structure

Table 2 -
Goudie (2006)n of rock-type, strength and fluvial erodibilityRock strength classification and descriptions modified from Table2inGoudie (2006) 4 ResultsIn plan view(Fig 1a), the drainage divide follows a sinuous form, with lengths mostly configured to ENE-WSW and WNW-ESE orientations.The divide tends to occupy central positions in relation to the topography of the range apart from a notable southerly segment in the middle of the study area.