Bedload-Bedrock Contrasts Form Enigmatic Low-Relief Surfaces of the Pyrenees

Low-relief, high-elevation surfaces

side of the range are well correlated up to approximately 1,800 m, with small remnants continuing up to 2,900 m (Babault et al., 2005;Bosch et al., 2016). On the northern side, a similar concentration of surfaces are correlated up to around 800 m (Bernard et al., 2019). Interpolating a continuous surface through LRHES suggests that they potentially originate from a gently undulating landscape (Bosch et al., 2016).
Geomorphic and geodynamic arguments have been put forward to refute the two leading scenarios. For the first scenario, in which peneplanation occurred around 800 m and drained directly to sea-level, followed by a subsequent increase in rock uplift, there is limited evidence for enhanced crustal thickening or removal of mantle lithosphere below the central Pyrenees (Bosch et al., 2016). For the second scenario, in which the surface formed at high elevation in response to continued erosion and sediment aggradation along the Pyrenean piedmonts, it is unclear whether such a large amount of aggraded sediment could be efficiently evacuated without record (Gunnell & Calvet, 2006). Recent work investigating sediment transport from a flexural isostatically compensated orogen to a coupled foreland basin suggests that some sediment drape and aggradation is expected as seen on the northern Pyrenees, but not to the highest elevations of the observed surfaces (Bernard et al., 2019). In addition, the sediment load associated with this aggradation would lead to significant flexural deflection of lithospheric plates, yet the overall load of the central Pyrenean topography has remained relatively constant since the Oligocene (Curry et al., 2019).
Here we explore an alternative process which demonstrates how evolution of exposed lithology disrupts drainage networks and creates LRHES. By accounting for contrasts between bedload and bedrock erodibility, we show that even where bedrock strength contrasts are very large, channel steepness values may be relatively uniform across lithological units, thereby obscuring the apparent importance of lithology.

Rock Strength Estimates From the Pyrenees
In the case of the Pyrenees, the importance of lithology shaping topography is evidenced by the main drainage divide tracking many of the exposed crystalline basement massifs (Figure 1b; Bernard et al., 2019). Metrics derived to evaluate the erodibility of Pyrenean landscape show a weak correlation with lithology and metric values overlap across different lithologies (Bernard et al., 2019). One metric is based on a model of fluvial elevation through time, dz/dt, responding to rock uplift rate, u, and erosion rate, e. In the stream power model, e is the  (Lehner et al., 2008). This zone is referred to as the Axial Zone. High elevations are found within the Axial Zone however, there is no clear correlation with lithology. (b) The tortuous path of the main drainage divides appear to coincide with the locations of the plutonic rocks as quantified by Bernard et al. (2019). (c) An example of a debated low-relief surface, the Plan de Beret. This surface is at an elevation of ∼1,850 m and steep topography is identified below and above the surface. Image from Google Earth. (d) Differential erosion is observed across the Pyrenees as highlighted by this active incision into a relatively low-relief surface near to Nerin. Image from Google Earth. product of the upstream drainage area, A and the local channel slope, S, raised to the powers of m and n, respectively (Howard, 1994) and erodibility, K, which encompasses bedrock strength, bedload, hydraulic parameters and climate. Therefore, At close to steady-state, as may be expected for an inactive mountain range in which rock uplift is driven by isostatic compensation to erosion like the Pyrenees, dz/dt ∼ 0, and u = KA m S n . This is commonly rearranged in terms of a normalized channel steepness index k sn = A m/n S, where m/n = 0.3-0.8 (Mudd et al., 2018), which provides an estimate of the spatial variability in u or K (Kirby & Whipple, 2012). A limitation with this approach is that maps of k sn can be very noisy because calculating slopes from DEMs accentuates noise. Therefore, methods are required to smooth and average the data so that coherent signals can be extracted from the noise. Here we use the approach of Fox (2019) to extract robust values of k sn (Figure 2a) as this approach is less sensitive to noise (Smith et al., 2022). k sn values are similar across the massifs and the surrounding rock ( Figure 2b). This is unexpected considering the correlation between k sn and lithology in the similarly post-orogenic southern Appalachian Mountains (Gallen, 2018). We also show that normalized channel steepness values calculated on a pixel-by-pixel basis show no clear correlation with lithology. Furthermore, the uplift rate values are expected to be similar across the Pyrenees and this suggests that the local erodibility, K, is varying. The relationship between compressive strength and tensile strength is roughly linear and therefore, despite the inherent limitations of the Schmidt Hammer, we assume that the tensile strength difference between the plutonic rocks and surrounding country rocks varies by a factor of two. Therefore, there is a clear discrepancy between the expected variability in the erodibility if compressive strength is a suitable proxy for erodibility.
The importance of bedload in driving river incision has been described mechanistically and measured in laboratories (Sklar and Dietrich, 2001). These experiments showed that the grainsize and sediment flux control erosion rate and this has been supported by field data (Brocard et al., 2016;Callahan et al., 2019;Finnegan et al., 2017;Johnson et al., 2009;Shobe et al., 2018). However, there is no evidence that grainsize varies with lithology within sedimentary rocks deposited at the time of Pyrenean Orogenesis (Michael et al., 2013).
In summary, the four key observations that motivate our model are: (a) patchy low-relief, high-elevation topography with evidence for transient incision; (b) drainage divides are tortuous and generally coincident with harder plutonic rocks; (c) contrasts in Schmidt Hammer measurements between the two dominant rock types; and (d) k sn variations are comparable between the two dominant lithologies. Our model is based on allowing K to evolve through time as a function of the contrast between the bedrock and bedload erodibility values. This means that the effective erodibility is very variable within a lithological unit, obscuring k sn -lithology variations. In turn, drainages networks are disrupted due to evolving short-wavelength changes in local erosion rates, producing tortuous divides that tend to coincide with the lower erodibility plutons.

A Simple Model Accounting for Bedrock-Bedload Contrasts
The expected erosion due to the impact between a particle being transported by a river and the underlying bedrock can be written as a function of the bedrock and particle strength (Sklar and Dietrich, 2004). In particular, the resistance of bedrock to erosion by impacting particles, , can be written as where is the tensile strength of the bedrock, Y is the Young's Modulus of elasticity of the bedrock and is a dimensionless coefficient that depends on the material properties of the impacting particle (Engle, 1978). For our illustrative model, in which we are trying to explore the impacts on how changes in outcrop through time may change downstream erodibility, we make the simplifying assumption that the channel erodibility is inversely proportional to this resistance. So, the erodibility is given by 10.1029/2022GL101995 4 of 10 where K c is a constant of proportionality and includes the Young's modulus, which is assumed to be constant.
Based on the Schmidt Hammer by Bernard et al. (2019), the tensile strength difference between the plutonic rocks and surrounding country rocks varies by a factor of two. If erodibility is proportional to 1/ 2 (Equation 2), the  (2019). This method regresses a variable channel steepness map through the relationship between χ and elevation keeping the baselevel for each catchment fixed. This result is produced using a pixel size of 2 km × 2 km. Half of all fluvial nodes with an upstream drainage area greater than 5 km 2 were randomly selected for the simplified topographical data set inversion resulting in a total of 77760 data points. A value of m = 0.45 and a smoothing parameter of 5 were used. (b) The distribution of normalized channel steepness values for pixels in the simplified data set. At each node in the simplified data set a channel steepness value from the inversion is extracted. erodibility of the plutonic bedrock is approximately four times less than the erodibility of the country rock. We term this contribution to the erodibility and the plutonic erodibility is ( ) and the country rock erodibility is ( ) . In addition, based on the similarity between channel steepness values across the granitic and county rocks, we can say that the erodibility for the plutonic rocks with plutonic bedload is equal to the country rocks with country rock bedload: which can be rearranged to give, The simplest solution to this is that ( ) is equal to ( ) and that ( ) is equal to ( ) . This suggests that the material properties of the bedload can be taken as the erodibility of the bedrock. More generally we can say that the local erodibility K is equal to where K c is simply a scaling parameter. The value of is simply taken from the average erodibility of the bedrock upstream of a specific location, as shown in Figure 3. Similarly, K r is the erodibility of the local bedrock. In this way, we do not actually track the bedload but assume that it quickly equilibrates to represent the upstream bedrock lithology.
In order to simulate the progressive outcropping of plutonic rocks, oblate spheres of low bedrock erodibility are exhumed toward the surface (Scherler & Schwanghart, 2020), driving river network reorganization (Bernard et al., 2021). The size of our landscape evolution model is similar to the overall size of the Pyrenees. The rock uplift rate is set to 1 km/Ma as a compromise between the pulses of rapid exhumation during the Oligocene and the low exhumation rates since (Fitzgerald et al., 1999;Gibson et al., 2007;Gunnell et al., 2009). The constant erodibility parameter ( ) is chosen to result in approximately 2 km of relief as observed in the Pyrenees today. The number of plutons is chosen to result in a modern-day distribution similar to what is observed in the Pyrenees, and these are located randomly in time and space. The value of n is set equal to one to highlight the first order influence of accounting for bedload composition in drainage network reorganization during landscape evolution. We acknowledge this is a simplification but it is sufficient to highlight our basic concept. In particular, Perne et al. (2017) show that the value of n is not equal to one, the slope of the contact between lithologies can control local incision using elegant numerical models. However, at the scale of our model, all lithological boundaries are assumed to be vertical. The area exponent, m, is equal to one half, and we include a diffusion term in our landscape evolution model to account for hillslope processes resulting in a mass balance equation: where dz/dl is the local slope in the direction of maximum slope and ∇ 2 is two-dimensional curvature. The erodibility, K, is given by Equation 7 and is the diffusivity and is equal to 1 × 10 −2 m 2 /a. This model is solved across all parts of the landscape evolution model. We discretize space into pixels of size 1 km × 1 km, the model runs for 5 Ma with time steps of 0.01 Ma. We use the Fastscape algorithm to solve for elevations through time (Braun & Willett, 2013).

Results
The difference between the erodibility of the plutons (1.25 1/Ma) and country rock (5 1/Ma) is large (Figure 3b), and local variations in K can be even larger (Figure 3c) but a clear correlation between k sn and bedrock lithology is not observed. Evolving differences in the erodibility lead to local changes in erosion rate (Figure 4). For example, when a point on the river network within the country rock has plutonic rock upstream, the bedload is made up of the hard plutonic rock and the bedrock is made up of the softer country rock. This produces an effective erodibility value of up to 20 1/Ma which produces high local erosion rates. These erosion rate variations lead to elevation changes and drainage network reorganization. When the disruption of drainage networks results in reductions in upstream area, fluvial erosion decreases, yet hillslope processes continue to operate, thereby reducing local relief (Yang et al., 2015). The surfaces, highlighted by the areas of low erosion rate (Figure 4), are found in both the plutonic rocks and the country rock. Once these surfaces have formed, they are then incised into as has been observed in the Pyrenees (Uzel et al., 2019;Figure 1d). In the Pyrenees, the calculated drainage divides tend to follow plutons (Bernard et al., 2019), in our models ∼1% of the outcropping country rock compared to ∼5% of the plutonic rock has zero upstream drainage area. This shows that the plutonic rocks are more likely to form drainage divides without having a k sn signature of a much harder rock. Importantly, when bedrock differences in erodibility are used as the effective erodibility and bedload-bedrock The erodibility values used in the stream power model is a function of the local contrast between the bedload and the bedrock erodibility. At some locations low erodibility bedload is transported over country rock, resulting in high erodibility values. At other locations the opposite is true. The inset shows that the normalized channel steepness values for the two lithologies are similar and this is to be expected given the similarity in erodibility values.
contrasts are not accounted for, k sn values strongly depend on lithology. Furthermore, spatial variations in erodibility are reduced from varying between >1 and ∼20 1/Ma in Figure 3c to 1.25 to 5 1/Ma in Figure 5c. This reduction in the variability of erodibility reduces the spatial variability in erosion rates. Ultimately, spatial variability in erosion rates across drainage divides disrupts the drainage network and thus when erodibility only varies as a function of bedrock lithology, drainage divides are not so dynamic (Figures 5 and 6). This reduction in divide dynamism, shown by reduced low-erosion rate locations, reduces the chance of LRHES formation.

Discussion and Conclusions
Our simple model reproduces several of the key features of the topography of the high Pyrenees and does not require complex changes in boundary conditions (Figures 3 and 4). Instead, we rely on superficial changes in outcropping lithology and on steady tectonic forcing. This relatively steady tectonic forcing is supported by post-orogenic flexural modeling (Curry et al., 2019). Overall, drainage divides are most likely to be found within the low erodibility plutons in our model, as is the case with the massifs of the Pyrenees. Furthermore, bedload contrasts obscure trends between lithology and normalized channel steepness, despite large differences in relative strength between lithologies. The degree to which the locations of the plutons control the positions of the drainage divides depends on the time that the plutons appear at the surface and the erodibility contrasts imposed in the model. The size of the plutons, the background rock uplift rate and the erodibility contrast between the plutons and the country rocks will all influence these outcomes, but the simple conclusions of our analysis are robust across a suite of model configurations tested.
This new mechanism for the formation of LRHES is driven by lithological contrasts between bedload and bedrock erodibilities and the resulting drainage network reorganization. Traditional mechanisms to create these surfaces have required large scale geodynamic events such as peneplanation close to sea level followed by rock uplift (Gunnell et al., 2008) or km-scale changes in fluvial baselevel (Babault et al., 2005). In other locations, recent research has deviated from these concepts and suggests continued tectonic evolution can drive drainage network reorganization to produce LRHES but these have still been based on spatial changes in rock uplift (Yuan et al., 2022) or changes in tectonic shortening (Yang et al., 2015). In contrast, our mechanism is based on steady tectonic processes with evolving lithology, creating drainage divides pinned to plutonic rocks but no clear lithological signature in k sn values. More generally, our two-dimensional landscape evolution model highlights how bedload composition evolution sets incision rates and, ultimately, drives drainage network reorganization. It is therefore important that resulting k sn variations from processes related to changing upstream areas are not mistaken for signatures of lithology or tectonic control on river incision. In all panels, there are areas of low erosion rate, both large areas forming low-relief surfaces and small areas that are preserved from previous time steps. From 4.5 to 5 Myr, the small pluton in the north is exhumed and blocks off a large area. Here, fluvial erosion rates drop while hillslope diffusion continues. In addition, at 5 Myr two plutons in the south east have joined and the smaller of the two plutons is dragging the divide to the south, as shown by asymmetric erosion rates across the divide immediately to the north.

Figure 5.
If the effective erodibility is simply a function of the bedrock lithology, the lithology has a clear impact on the channel steepness. This is shown in the insert in (c). The layout of the figure is the same as Figure 3. The effective erodibility is less variable compared to Figure 3c and varies between 1.25 and 5 1/Ma. This change in erodibility leads to systematic patterns in erosion rate between the two lithologies producing distinct channel steepness signatures. The variability in the channel steepness values observed in the insert is the result of the evolution of outcropping lithology and drainage network reorganization.

Figure 6.
If the effective erodibility is simply a function of the bedrock lithology the overall topography is clearly a function of the lithology, with plutonic units forming the highest parts of the landscape. In this way, the drainage network is also very dependent on the lithology. In addition, the drainage network is less dynamic and the mechanism of lithological variations forming high elevation surfaces that cut through different lithological units is less effective due to the reduction in the range of effective erodibility. The layout is the same as Figure 4.