Investigation of stochastic-threshold incision models across a climatic and morphological gradient

. Long-term landscape evolution is controlled by tectonic and climatic forcing acting through surface processes. Rivers are the main drivers of continental denudation because they set the base level of most hillslopes. The mechanisms of ﬂuvial incision are thus a key focus in geomorphological research and require accurate representation and models. River incision is often modeled with a stream power model (SPM) based on the along-stream evolution of drainage area and channel elevation gradient but can also incorporate more complex processes such as threshold effects and statistical discharge distributions, which are fundamental features of river dynamics. Despite their importance in quantitative geomorphology, such model formulations have been confronted with ﬁeld data only in a limited number of cases. Here we investigate the behavior of stochastic-threshold incision models across the southeastern margin of the French Massif Central, which is characterized by signiﬁcant relief and the regular occurrence of high-discharge events. Our study is based on a new dataset combining measurements of discharge variability from gauging stations, denudation rates from 34 basins from 10 Be cosmogenic radionuclide (CRN) concentration measurements in river sediments, morphometric analysis of river long proﬁles, and ﬁeld observations. This new dataset is used for a systematic investigation of various formulations of the SPM and to discuss the importance of incision thresholds. Denudation rates across the SE margin of the Massif Central are in the 20–120 mm kyr − 1 (equivalent to mm/ka in the ﬁgures) range, and they positively correlate with slope and precipitation. However, the relationship with the steepness index is complex and supports the importance of taking into account spatial variations in parameters ( D 50 , discharge variability k , runoff) controlling the SPM. Overall, the range of denudation rate across the margin can mainly be explained using a simple version of the SPM accounting for spatially heterogeneous runoff. More complex formulations including stochastic discharge and incision thresholds yield poorer performances unless the spatial variations in bedload characteristics controlling incision thresholds are taken into account. Our results highlight the importance of the hypotheses used for such a threshold in SPM application to ﬁeld studies and notably the impact of actual constraints on bedload size.


Introduction
Investigating landscape evolution and its responses to tectonic or climatic perturbations requires accurate models of surface processes such as hillslope erosion and fluvial incision (Dietrich et al., 2003). Rivers are key agents for the transmission of external forcing through landscapes and, by setting the base level of most hillslopes, the main drivers of continental denudation. For these reasons the mechanisms of fluvial incision have been a key focus of geomorphological research over the last 3 decades (e.g., Howard et al., 1994;Whipple et al., 2000).
A powerful framework has been derived for the analysis of fluvial processes along bedrock rivers by postulating that incision rate scales with either shear stress or unit stream power per riverbed area (e.g., Howard, 1994;Whipple and Tucker, 1999;Sklar and Dietrich, 2006;Turowski et al., 2007;Lague, 2014). The main appeal of this stream power model (SPM) is that the incision rate can be expressed using a limited number of reasonable assumptions as a simple function of drainage area and channel elevation gradient, which are variables easily extracted from digital topographic data. Such a simplification is the foundation for the widely used analysis of channel steepness, which is now a standard approach in the investigation of tectonically active landscapes, allowing for deciphering relative rock uplift distributions (Wobus et al., 2006;Kirby and Whipple, 2012;Whittaker, 2012). It has also been used to model fluvial incision at continental scales (e.g., Roberts and White, 2010). However, the tremendous success of these applications of the stream power model should not hide the fact that some of its key aspects remain highly debated in terms of both the parameters it incorporates and the physical processes it is supposed to represent (e.g., Lague, 2014;Gasparini and Brandon, 2011;Harel et al., 2016). For example, the way rock resistance to erosion can be accounted for is still poorly understood (e.g., Sklar and Dietrich, 2001), as is its applicability to fluvial systems wherein sediment transport and deposition represent an important modulation of incision (e.g., Whipple and Tucker, 2002).
The existence of a threshold in discharge that needs to be overcome for incision or sediment transport to occur has long been recognized in fluvial geomorphology. Taking into account such a threshold in relation to a statistical description of discharge is a long-standing issue for the modeling of long-term fluvial incision and for our understanding of the contribution of extreme events to landscape evolution (e.g., Tucker and Bras, 2000;Snyder et al., 2003;Lague et al., 2005;Molnar et al., 2006;Deal et al., 2018). The theoretical implications of incorporating incision thresholds and discharge distributions into stream-power-based incision models have been explored systematically and highlight the emergence of complex and nonlinear behavior between channel properties and incision rates.
Recently, several studies have built upon this theoretical framework and explored its applications to real landscapes, mainly by bringing together erosion rates derived from cosmogenic radionuclides (CRNs), river profile morphological analysis, and constraints on climatic or hydrological variability. They notably highlighted the importance of taking into account incision thresholds and discharge variability to understand the nonlinear relationship between erosion rates and channel steepness. DiBiase and Whipple (2011) provided the first of such investigations in the San Gabriel Mountains in California and demonstrated the suitability of river incision models combining an incision threshold and a power-law distribution of extreme discharge events to account for observations of river profile steepness and erosion rates. More generally, they explored the behavior of this relationship across a wide range of climatic configurations to highlight limits in the response of channel incision to an increase in runoff. Based on a denudation dataset spanning very contrasting climates along the central Chilean Andes, Carretier et al. (2013) demonstrated the role of extreme discharge events in driving most of the erosion for arid settings. Scherler et al. (2017) developed an analysis at the scale of the whole Himalayan arc and eastern Tibet in order to explore the underlying reasons for long-wavelength variations in the relationship between channel steepness and erosion rates. They notably explored the influence of a spatially variable incision threshold and assessed the importance of climatic controls such as characteristics of monsoon precipitation. More recently, Campforts et al. (2020) highlighted the importance of taking into account lithological and climatic spatial differences when using this type of model to analyze the variance in CRN-derived denudation rate datasets.
These studies were critical tests of our ability to reconcile topographic and geochronological observations using a physically based framework for river incision. However, many of these earlier investigations relied on pre-existing CRN datasets that were acquired with sometimes independent sampling strategies. It implies datasets combining catchment areas varying over 2 orders of magnitude and with notably large basins (> 100 km 2 ), which cover surfaces displaying important variations in bedrock geology and climate. Additionally, due to the spatial extent of the studied regions and the lack of appropriate gauging station data, some of these studies also relied in part on remote sensing observations for the characterization of hydrological variability instead of direct measurements. Lastly, in some cases it was not possible to constrain first-order parameters with field observations, such as bedload size, which is the most important factor controlling the incision threshold.
The main objective of our study is to investigate stochastic-threshold incision models in a context in which robust, short-wavelength constraints can be obtained for all the parameters involved. In particular, we want to explore the importance of spatial variations of environmental parameters such as incision threshold and discharge variability in order to better understand their control on long-term landscape evolution. In contrast with most previous studies which have addressed such questions, we design a dedicated sampling strategy focusing on the southeastern margin of the Massif Central (France), a region with well-defined shortwavelength climatic and hydrological gradients and limited recent tectonic activity. For that purpose, we acquired a new consistent CRN dataset for catchment-wide denudation to compare measured landscape evolution rates with predictions of various river incision models. We analyzed river profile and basin morphometric parameters to evaluate their relationships with denudation rates, compiled hydrological data from gauging stations in order to characterize discharge variability, and combined new and existing bedload size data to constrain the incision threshold.
Based on the premise that river networks are the main drivers of landscape evolution, the relationship between channel steepness and erosion rates has been intensively scrutinized, in particular to assess the implications of its degree of nonlinearity. At first order, our new dataset displays a complex structure between these two observations, which cannot be adequately fitted by a simple scaling relationship, and highlights the need for a careful consideration of spatial variations in the lithological, climatic, and hydrological parameters. We develop this argument by first recalling the theoretical background behind stochastic incision models and outlining the main characteristics of our study area, with an emphasis on the specific characteristics which make it a particularly appropriate setting for our investigation of such models. Then, we describe the methods used and the dedicated sampling strategy we designed, as well as the results obtained. On the basis of these results, we explore the consistency of model predictions with the observed denudation rates. Lastly, we discuss the regional distribution of these denudation rates, their co-variations with various environmental factors, the adequacy of various modeling approaches, and the implications for stochastic incision models.

Theoretical background
The theoretical background for the incorporation of discharge probability distributions into river incision models has been extensively described in previous studies (Tucker and Bras, 2000;Snyder et al., 2003;Lague et al., 2005;DiBiase and Whipple, 2011;Scherler et al., 2017;Campforts et al., 2020), and its main features will be presented briefly here. Instantaneous detachment-limited fluvial incision I can be expressed as a function of shear stress acting on the riverbed τ , where τ c is a threshold shear stress beyond which incision occurs, k e is an erodibility coefficient depending on the properties of the bedrock, and a is an exponent dependent on the nature of incision processes. Shear stress τ can be expressed as a function of discharge Q, slope S, and channel width W as with k t , α, and β as parameters depending on the formulation used to express flow resistance: for example, α = 3/5 and β = 7/10 when using Manning's frictional relationship or α = β = 2/3 for a Darcy-Weisbach relation (Howard, 1994). The critical shear stress can be expressed as where D 50 is median bedload grain size, τ * c is the critical Shields stress, and ρ s and ρ w are the density of sediment and water, respectively. The variation of channel width W downstream and with instantaneous discharge Q can be computed according to with ω s an empirical parameter and W b the bankfull width, scaling with the mean annual discharge Q as where k w and ω b are another couple of empirical parameters, which can be determined from field measurements or remote sensing imagery (Kirby and Ouimet, 2011;DiBiase and Whipple, 2011;Fisher et al., 2012). Mean annual discharge Q can be calculated from mean runoff R and drainage area A as Q = RA c , with c an exponent that we will consider equal to 1 in the following. Equations (1) to (5) can then be combined into an expression for instantaneous incision I , This equation can be simplified by defining an area exponent m = aα(1−ω b ), a slope exponent n = aβ, an erodibility term K = k e k t k −aα w R m , γ = aα(1−ω s ), and a threshold term = k e τ a c , This relationship highlights the importance of the parameter γ in controlling the sensitivity of incision to discharge and its variability (Lague et al., 2005). Introducing the steepness index k s = A m/n S (Kirby and Whipple, 2012), Eq. (7) can also be written as Assuming that the threshold term is negligible and that discharge Q is constant and equal to mean discharge Q, Eq. (7) 476 C. Desormeaux et al.: Investigation of stochastic-threshold incision models simplifies to the classical stream power form, which is widely used in quantitative geomorphology studies: Long-term incision E can be calculated by weighting instantaneous incision (Eqs. 7 or 8) with discharge probability pdf(Q) and integrating over Q, where Q c is the critical discharge needed to overcome the incision threshold and Q m is the maximum discharge. Q c is obtained by setting I = 0 in Eq. (7): In incision models taking into account discharge variability, the probability distribution of discharges is often modeled using an inverse gamma law, where Q * = Q/Q is normalized discharge, k is the variability parameter, and is the gamma function (Lague et al., 2005). Low values of k correspond to a situation in which the discharge variability is high. The inverse gamma distribution combines exponential and power-law parts describing low-and high-discharge regimes, respectively, and accounting for the low probabilities of events at both ends of the spectrum. In this study we will use Eq. (10) to model erosion rates using field constraints on its parameters and compare its predictions with cosmogenic-nuclide-derived denudation rates.

Field setting
Our region of interest is located along the southeastern margin of the French Massif Central, directly west of the Rhône river valley (Fig. 1a).
The margin presents an asymmetrical topography with deeply incised valleys draining into the Rhône on its eastern flank (Fig. S4), while farther west, the interior of the Massif Central is characterized by lower relief and plateau areas (Fig. 2).
The limit between these two domains corresponds to the divide between the Rhône and Loire or Garonne catchments, and on the Rhône side of the divide we specifically investigate the Cévennes and Ardèche mountains with maximum elevation of ∼ 1600 and 1750 m, respectively.
Orographic precipitation associated with moisture from the Mediterranean is focused on the topographic margin, with mean annual precipitation (MAP) ∼ 1500 mm, and then a slight decrease toward MAP ∼ 1000 mm inside the Massif Central (Figs. 1a and 2). This orographic precipitation regime is also associated with very intense rainfall events, known as Cévenol events, triggering flash floods along the rivers draining the southeastern margin, usually in the fall. For example, during the 8-9 September 2002 event ∼ 500 mm precipitation was recorded in 9 h at one rain gauge inside the Gardon catchment (Le Lay and Saulnier, 2007).
The Massif Central is part of the Hercynian orogenic system, which developed during the late Paleozoic. Due to their connection to the Rhône valley the eastern drainage systems were probably affected by the Messinian salinity crisis (Mocochain et al., 2009;Tassy et al., 2013). The Cenozoic uplift history of the Massif Central and its mechanisms are still poorly documented, but it is suspected that most of the observed relief, with respect to the Rhône valley, is associated with distinct late Miocene to Pliocene uplift events (Olivetti et al., 2016;Malcles et al., 2020). The present tectonic activity of the Massif Central is limited, with a recent uplift rate < 100 mm kyr −1 . A few historical and instrumental earthquakes with magnitude < 5 have been reported along the southeastern margin (Mazzotti et al., 2020;Ritz et al., 2020). Bedrock lithology is dominated by crystalline Hercynian basement (granites, gneisses, and micashists), Mesozoic-Cenozoic sedimentary series, and late Cenozoic volcanic rocks (Fig. S1 in the Supplement). In the following we focus on the quartz-bearing lithologies when collecting samples for 10 Be analysis.
In summary, the SE Massif Central presents a clear topographic margin with contrasting domains on both sides of the drainage divide, associated with strong climatic and hydrological gradients. In our study we take advantage of this particular configuration, as well as low tectonic activity and strong climatic variability, to explore stochastic-threshold incision models.

Methods and data
We present the various datasets which were acquired or used for our study. For the sake of clarity, for each type of data we successively present the methodology and corresponding results.

Methods
We extracted daily discharge data for 326 gauging stations from the Banque HYDRO (http://hydro.eaufrance.fr, last access: 31 May 2022) database over the SE margin of the Massif Central and surrounding areas, with records spanning at least 20 years. The time series were first processed to identify anomalies such as missing values, zero values, or sequences of days with identical data. For each station starting and , with the location of the sampled basins and their denudation rate (red numbers). Black numbers identify the different samples. The red lines are the main regional drainage divides between the majors catchments: Loire (Lo), Rhône (Rh), and Garonne (Gr). For basins in the Rhône catchment, we distinguish two areas corresponding to the Cévennes (Rh-C) and Ardèche (Rh-A) mountains (Table S3). (b) Distribution of precipitation over the studied area. The black rectangle indicates the extent of (a). Colored circles are the gauging stations where the discharge variability coefficient k was estimated, with dark contours indicating selected stations with discharge distributions presented in Fig. 3. Dark contours correspond to a thin plate spline surface fitted to the variability data. aa' indicates the position and lateral extent of the swath profile presented in Fig. 2. ending dates were manually defined in order to extract the longest continuous series of valid data. The discharge data were normalized by the mean discharge and empirical cumulative distributions when calculated. Our focus here is on the high-discharge part of the distribution and the contribution of events associated with incision, and we fitted a power law to that part of the distribution of the form pdf(Q) ∝ Q −α . The value of exponent α is an indicator of the importance of highdischarge events in the flow regime (Molnar et al., 2006) and is directly related to the variability parameter of the inverse gamma distribution (Eq. 12) by α = k + 2.
The power law was fitted using the approach of Clauset et al. (2007) to determine both the exponent α and the minimal normalized discharge value corresponding to the lower bound of the power-law behavior. The identification of the lower limit x min of the power-law behavior was done by selecting the value yielding a probability distribution of the empirical normalized discharge data above x min as close as possible to the best-fit power-law model. In addition to a longwavelength pattern associated with the global regional climatic setting, we expect to observe short-wavelength variations due to diversity in catchment size, bedrock geology, and vegetation as well as uncertainty resulting from the limited length of the records, differences in monitoring practices by operators, and biases in the measurement of extreme events. We are interested in the first-order regional variation of climatic parameters; in order to obtain a long-wavelength approximation of the spatial distribution of k over our study area we fitted a thin plate spline surface to the station measurements of k.

Results
The calculated variability parameter k (Eq. 12) ranges from 0.1 to 3.7, with lower values similar to contexts characterized by pronounced discharge variability and a significant occurrence of extreme events (Lague et al., 2005;Molnar et al., 2006;DiBiase and Whipple, 2011). We observe a clear gradient in the spatial distribution of this variability parameter (Figs. 1b, 2, and 3).
Low k values are dominant along the topographic margin of the SE Massif Central in the Cévennes and Ardèche areas, but also farther south and north along the Rhône valley. It corresponds to the areas affected by intense precipitation events of Cévenol and Mediterranean type. Conversely, when mov- ing inside the massif, the measured k increases, with values in the 2 to 3.7 range indicative of a less variable hydrological regime. Furthermore, the thin plate spline surface correctly reproduces the global spatial variability gradient across the margin (Fig. 2), and we use it to compute basin-averaged k and incorporate it into the stochastic incision models.
Paleofloods and historical floods have been documented in our study area, in particular along the Ardèche river (Sheffer et al., 2003;Naulet et al., 2005) and the Gardon river (Dezileau et al., 2014), allowing the assumption of a high-discharge-variability regime since at least the middle Holocene. The absolute values of precipitation and discharge variability have probably changed through time and been modulated by the succession of glacial cycles, but, as the present-day relief is associated with late Miocene to Pliocene uplift events (Olivetti et al., 2016(Olivetti et al., , 2020Malcles et al., 2020), it can be assumed that the orographic pattern and associated relative climatic contrasts across the margin are persistent throughout the Quaternary.
While a number of previous geomorphological studies have investigated the impact of variations in parameters quantifying discharge variability over large areas such as the Himalayas (Scherler et al., 2017) and the continental US (Molnar et al., 2006), our study is the first to tackle this issue over a such a well-defined short-wavelength (< 100 km) gradient.

Methods
We only have discharge gauging station data at a limited number of locations, which are not coincident with our sampling sites. We calibrated a relationship between basinaveraged precipitation and runoff at these instrumented locations in order to be able to compute runoff at the sites where we obtained 10 Be denudation rates. We extracted catchment boundaries for Banque HYDRO stations in our study area and computed the catchment-averaged mean annual precipitation (MAP) using a 250 m resolution raster (Joly et al., 2010), as well as the runoff from discharge data for these stations. We restricted our analysis to drainage areas ranging from 10 to 400 km 2 .

Results
We observe a positive correlation between basin-averaged MAP and runoff measured at gauging stations (Fig. 4).
There is no clear distinction in behavior between the interior of the massif and its southeastern margin, and catchments draining over carbonate lithologies do no display a different trend than the rest of the dataset. We use the fitted linear relationship between MAP and R in order to compute an estimate for runoff R for the catchments where we obtained 10 Be denudation rates (Table S3).

Channel width 4.3.1 Methods
We calibrated the scaling relationship between discharge and channel width (Eq. 5) with measurements at various locations in our region of interest. We use IGN orthophotos (BD ORTHO IGN -50 cm resolution) in low-vegetation areas or field measurements with a laser range finder. We combined our data with measurements from the Carhyce database (Gob et al., 2014) (http://194.57.254.11/IED/, last access: 1 May 2021). Discharge values at the measurement site were calculated from the runoff-precipitation calibration relationship (Fig. 4). We restricted our analysis to drainage areas < 400 km 2 and discharge > 0.1 m 3 s −1 .

Results
We observe a positive correlation between channel width and discharge (Fig. 5a), and, for a given discharge, channel widths are larger along the margin (Cévennes and Ardèche) than west of the divide (red symbols).  We fit power laws to the whole dataset and to subsets corresponding to both sides of the divide in order to calibrate the relationship described by Eq. (5). When considering the whole dataset we obtain a value for the exponent ω b = 0.44, which is within the range of reported previous observations (Campforts et al., 2020) and close to the ω b = 0.55 we impose in order to keep Eq. (8) dimensionally consistent when using m/n = 0.45. In this case we obtain a value of k w = 5.56 ± 0.37 m −0.65 s 0.55 . If we distinguish the regions east and west of the divide, the values of k w are 4.95 ± 0.59 and 6.85 ± 0.49, respectively. We consider this difference in k w to be significant over our area of interest, and we use distinct values for the corresponding regions in the following.

River long profile 4.4.1 Methods
We extracted the catchment contour and upstream river network for each site where 10 Be concentrations were measured in order to compute basin-scale topographic metrics such as normalized channel steepness k sn (Kirby and Whipple, 2012). We use a 25 m resolution digital elevation model (BD ALTI IGN) and a source area of 2.2 km 2 for channel initiation. The optimal concavity of the stream network for each basin was defined using the integral approach of Perron and Royden (2012). Normalized steepness indexes were  Fig. 4. Circles, triangles, and squares correspond to orthophotos measurements, field measurements, and data from the Carhyce database (Gob et al., 2014), respectively. The dashed black line is a power-law fit of the whole dataset, whereas blue and red dashed lines correspond to fits to the Ardèche-Cévennes (high relief, east of the divide) and Loire-Garonne (lower relief, west of the divide) subregions, respectively. Symbols are colored according to the four different subregions. (b) Median bedload size D 50 as a function of steepness index k sn . Triangles and squares correspond to our field measurements and the Carhyce database, respectively. Yellow squares are binned values (median and interquartile range), and the yellow line is the corresponding power-law fit (Eq. 13), with its 99 % confidence envelope. The dashed green line corresponds to the estimate of the parameters in Eq. (13) from the model presented in Fig. 7f. then computed along the river profiles that were filtered to remove short-wavelength noise and artifacts using a reference concavity value set at 0.45 (Fig. S3).

Results
The channel concavity ranges from 0.13 to 0.87 (Table S3) with a mean value of 0.52 ± 0.19, which is close to the usual reference concavity value used in channel profile analysis. The normalized steepness index k sn ranges from 22 to 97 m 0.9 . The highest values are observed in the Ardèche, Cévennes, and Garonne areas with averages of 63 ± 13, 47 ± 11, and 44 ± 12 m 0.9 , respectively, while the Loire watershed presents slightly lower values at 30 ± 7 m 0.9 .

Methods
We use bedload median size D 50 to estimate the critical shear stress τ c using Eq. (3). We performed bedload counts at some sites in the Cévennes and Ardèche areas and combined these results with D 50 measurements from the Carhyce database (Gob et al., 2014) (http://194.57.254.11/IED/, last access: 1 May 2021). In addition to the determination of a single τ c value derived from the average D 50 over our study area (DiBiase and Whipple, 2011), we also consider the possibility of spatial variations in τ c . Such variability can be parameterized using Eq. (3) and a power-law relationship between k sn for the upstream basin and the median bedload grain size D 50 (Attal et al., 2015;Scherler et al., 2017), such as

Results
Median bedload size ranges from 10 to 200 mm over our study area. The average D 50 is 80 mm (Fig. 5b), corresponding to τ c = 40 Pa, which is similar to the value used by DiBiase and Whipple (2011) in the San Gabriel Mountains. Coarse bedload was dominant in the river we surveyed along the margin, whereas we observed occurrences of gravel-or sand-dominated reaches for the low-relief parts of the landscape west of the divide. Fitting Eq. (13) to the dataset supports the existence of a power-law relationship between D 50 and k sn , with exponent ∼ 0.8. We do not observe a difference in scaling between the two sides of the divide. For the lowk sn regions, corresponding mostly to areas west of the divide, measured D 50 values are dominantly in the 30-60 mm range, whereas they are mostly > 50 mm for the eastern catchments draining toward the Rhône.

Methods
We measured 10 Be concentrations in quartz from river sands of basins draining the SE margin of the Massif Central, as well as the upper Loire and Garonne catchments. Sampled  8) and (10) for various values of the discharge variability parameter k. We set the critical shear stress τ c = 40 Pa and the erodibility coefficient k e = 8 × 10 −13 m 2.5 s 2 kg −1.5 . We use a regionally averaged value for runoff of R = 650 mm a −1 .
basins were primarily selected according to their size (21 to 95 km 2 with the exception of basin CDX-30) and the observation of regular concave river long profiles, without major knickpoints (Fig. S3). We also restricted our sampling to basins draining uniform lithologies, focusing on the Paleozoic basement units, constituted of granites and gneisses, as well as schists in the Cévennes area (Fig. S1). Quartz was extracted and purified from the samples by standard physical and chemical procedures. A mass of ∼ 20 g of quartz was digested with hydrofluoric acid (HF) along with a 9 Be carrier solution, and Be was extracted with ions exchange columns. 10 Be/ 9 Be ratios were finally measured at the French accelerator mass spectrometer (AMS) national facility, located at CEREGE in Aix-en-Provence. Steady-state denudation rates were computed using the online calculator described in Balco et al. (2008). A detailed description of the analytical procedure, calculation, and results is provided in the supplementary materials.

Results
Measured 10 Be concentrations in the 34 sampled catchments range from 51 to 345 × 10 3 at g −1 (Table S1 in the Supplement). The corresponding calculated denudation rates range from 24 to 126 mm kyr −1 (Fig. 1a), which correspond to integration times of 25 and 5 kyr, respectively. We observe clear regional differences in denudation rates between the main sampling areas. The upper Loire and Garonne catchments, in the interior of the Massif Central, display the lowest denudation rates with averages of 33 ± 4 and 37 ± 4.5 mm kyr −1 , respectively. Significantly faster denudation is observed along the margin of the massif, in the Ardèche mountains (58±7 mm kyr −1 ), and especially in the Cévennes area (95 ± 11 mm kyr −1 ). We note that basin CDX-14, with the highest denudation rate of the dataset, is located next to the main divide between the Garonne and Rhône watersheds, very close to the Cévennes area, and drains over schist bedrock. This range of denudation rates is consistent with previous observations in the Massif Central (Schaller et al., 2001(Schaller et al., , 2002Molliex et al., 2016;Olivetti et al., 2016). Denudation rates are positively correlated with basinaverage slope and relief (Figs. 6a and S2), but the dataset as a whole does not display any clear relationship with steepness index (Fig. 6b).
With the exception of CDX-CRN-30, all the sampled basins have areas between 20 and 100 km 2 , and we do not observe any dependence of denudation rates on catchment area over this range ( Fig. S2 and Table S3), suggesting that there is no first-order scale-dependent bias in terms of the processes contributing to denudation in our dataset (Niemi et al., 2005;Yanites et al., 2009).

Modeling approach
We model denudation rates for each basin where we have 10 Be data using various formulations of the stream power model. First, if we do not take into account discharge variability, denudation rate E can be derived from Eq. (9) as E = Kk n s = k e k t k −aα w R m k n s .
In this case we can either neglect the spatial variations of runoff, consider R to be constant, and set it to a regional av- erage value for all the considered basins (area-based stream power model: A-SPM) or take into account these spatial variations (runoff-based stream power model: R-SPM). In both cases, all the parameters are set to values from the literature or calculated for the areas we investigate (Table 1). The only unconstrained parameter is the erodibility coefficient k e . We account for differences in erodibility between the Cévennes schists and the granitic or gneissic basement encountered for all other areas by introducing a modulation factor f s for basins draining over the former. Thus, the A-SPM and R-SPM cases have two free parameters k e and f s .
We also consider stochastic-threshold stream power models (ST-SPMs), incorporating the probability distribution of stream discharge as formulated in Eq. (10). First, we set the critical shear stress τ c (Eq. 1) using the average observed median grain size over our study area using Eq. (3) (Fig. 5b). Second, similarly to Campforts et al. (2020) and Scherler et al. (2017), we consider τ c to be a free parameter in addition to k e and f s .
The last model formulation we consider is the case in which the critical shear stress τ c is not uniform over the study area, as in the previous case, but spatially variable (SVT-SPM) and defined using Eqs. (3) and (13), as explored by Scherler et al. (2017). Contrary to some previous studies (Scherler et al., 2017;Campforts et al., 2020) we have actual field constraints on D 50 and its relationship with k sn , and we first set the parameters k 50 and q of Eq. (13) according to these constraints (Fig. 5b). In a second run the models consider k 50 and q to be free parameters, similar to the approach of Scherler et al. (2017).
Following Campforts et al. (2020) we use various metrics to evaluate the quality of the fit between the n modeled (M) and observed (O with uncertainty σ ) denudation rates, such as the coefficient of determination R 2 for a linear fit between O and M or the χ 2 value, and we maximize Nash-Sutcliffe (NS) model efficiency for our determination of optimal parameters,

Results
We observe that R-SPM allows for a slightly better prediction of denudation rates than A-SPM ( Fig. 7a and b), with notable differences between the two models for lower denudation rates (< 50 mm kyr −1 ). In both cases the erodibility coefficient k e for the Cévennes schist bedrock is higher than for the crystalline basement lithologies of the other areas, which is consistent with estimates of relative differences in erodibility for these types of rocks (Campforts et al., 2020). Using an ST-SPM, with a unique threshold value calculated according to observed regionally averaged D 50 (Fig. 7c), induces a drastic degradation of the fit. While the model is able to reasonably predict the denudation rates for some of the high-relief catchments along the margin (Cévennes and Ardèche areas), it largely underpredicts them for the other areas inside the massif. The introduction of a variable threshold set by Eq. (13) (SVT-SPM) yields a better fit than the constant threshold case (Fig. 7d), notably by improving the prediction for some of lowest denudation catchments west of the divide. However, the quality of the fit is still inferior to the R-SPM case. When considering a free regionally constant threshold the optimization converges toward an almost negligible value, corresponding to a situation which is essentially similar to R-SPM (Fig. 7f). At last, if we consider the two parameters of the scaling described by Eq. (13) to be free parameters, the optimization yields a fit to the data which is equivalent to R-SPM. Noticeably, the retrieved scaling exponent q between D 50 and k sn of ∼ 0.8 is almost identical to the observed value (Fig. 5b), while the k 50 factor implies ∼ 2 times lower predicted D 50 values, but which are still within the lower range of observations.

Discussion
After the presentation of the data and main results of our study, we now discuss our denudation rate dataset in the context of the SE margin of the Massif Central, the relative efficiency of the various SPM formulations which can be used to explain the distribution of these denudation rates, and the comparison of our observations with similar studies in other regions. We note that the timescale of our study is relatively short and set by the integration times of denudation rates ranging from 5 to 25 kyr (Table S1).

Distribution of denudation rates
Here we discuss the spatial distribution of denudation rates across the SE margin of the Massif Central and their relationship with morphological parameters. As we described earlier, the range of denudation rates from our sampled catchments is consistent with previous observations in the Massif Central. Olivetti et al. (2016) report denudation rates ranging from 34 to 79 mm kyr −1 on 27 watersheds located directly north of our studied catchments in the Ardèche mountains, with similar lithologies and areas. As many of the catchments from Olivetti et al. (2016) contain knickpoints, the relationships between denudation rates and morphology is not directly comparable with the one observed with our dataset, for which we specifically selected basins according to the regularity of river long profiles. A larger-scale analysis of denudation in the Gulf of Lion drainage system was performed by Molliex et al. (2016), with four sampled watersheds along the southeastern margin of Massif Central, corresponding to the Eyrieux, Ardèche, Cèze, and Gard rivers near the Cévennes mountains, with average denudation rates of 79.6±15, 67.4±11.5, 78.0±16.2, and 70.0±11.9 mm kyr −1 , respectively. These basins are 1 order of magnitude larger than the ones we sampled and drain over heterogeneous lithologies with a significant contribution of the sedimentary cover in the low-relief area between the Massif Central margin and the Rhône river, which might explain the slightly lower denudation rates when compared to our samples from the Cévennes area. Denudation rates of 33 ± 7 and 68 ± 16 mm kyr −1 were obtained by Schaller et al. (2001) on the other side of the divide for the Allier and Loire catchments, respectively. Both catchments are larger (> 250 km 2 ) than the ones we sampled and include Cenozoic volcanics, but the denudation rate reported for the Allier headwaters from Schaller et al. (2001) is consistent with our observations from smaller catchments 11 and 12, while we have no equivalent in our dataset for the Loire headwater catchment where they report higher denudation rates. The current geography of the southeastern margin of the Massif Central and the observed contrasts in denudation rates suggest that there is some degree of westward migration of the divide as well as local river capture events, especially in the Cévennes and southern Ardèche areas. However, no large-scale disequilibrium such as knickpoints or relict surfaces associated with such captures was identified in the headwaters of our sampled basins. We also note that the measured concavities in our basins (Table S3) show a dispersion with relatively low values for some of the Loire and Garonne catchments (< 0.4) and higher values for some of the Cévennes and Ardèche catchments (up to 0.8 locally), although the averaged value is around 0.5, which is consistent with global compilations (Harel et al., 2016) and theoretical predictions of the stream power model with m = 0.5 and n = 1 when assuming a steady state over the study area (Lague, 2014). This observed dispersion in concavity may be due to the strong precipitation gradient over our study area (Fig. 2), but the correlation between denudation rates and concavity in some basins (Fig. S2) could also be interpreted as reflecting a transient response and retreat of the topographic margin. This deviation from a strict steady-state situation could explain some of the scatter in our dataset with respect to the predictions of the various modeling formulations we test.
Denudation rates are positively correlated with basinaverage slope (Fig. 6a), with a progressive linear increase in low-relief basins of the interior of the massif with slope ranging from 5 to 15 • and then a much more rapid increase beyond 20 • for the Cévennes and Ardèche basins along the margin. Such a relationship between denudation rates and slope is consistent with nonlinear models of hillslope evolution, whereby denudation rates increase very rapidly when hillslope angles become close to their stability threshold (Ouimet et al., 2009;DiBiase and Whipple, 2011;Carretier et al., 2013). Denudation rates are increasing nonlinearly with relief and also display a strong spatial difference between the margin and the interior of the massif (Fig. S2), with low denudation in the low-relief headwaters of the Loire and Garonne watersheds and significant scatter with higher rates for the Cévennes and Ardèche basins. We note that denudation rates are also positively correlated with runoff ( Fig. 8a). A similar correlation has been observed in other non-tectonically active contexts, such as Kaua'i island (Ferrier et al., 2013b, a), across a wide range of precipitation. As suggested by Ferrier et al. (2013b), lithological homogeneity and the absence of rock uplift gradient might be important factors in allowing the expression of such a relationship between denudation rates and precipitation. These two conditions are present in our study area; however, we note the existence of a pronounced orographic effect across the margin such that the correlation between runoff and denudation also partly reflects the underlying spatial distribution of precipitation with higher rainfall focused on the steep and high-relief areas (Figs. 1b, 8b and 2).
We do not observe a simple relationship between normalized channel steepness index and catchment denudation rates (Fig. 6b). Samples from the interior of the Massif Central display a large range of k sn up to ∼ 60 m 0.9 , while denudation rates remain low. In the Ardèche mountains a positive correlation exists between k sn and denudation rates, with a much shallower trend than for the Loire and Garonne catchments. Samples from the Cévennes area cluster outside the general trend, with high denudation rates for comparatively low k sn , which could in part be related to a lithological control on erodibility and the dominance of schist bedrock as opposed to the granitic and gneissic basement of the other sampled basins. Equations (8) and (10) predict a nonlinear relationship between denudation rate and steepness index k sn with a strong influence of the discharge variability parameter k (DiBiase and Whipple, 2011). Due to the absence of a clear unique trend in our dataset, no single theoretical relationship can explain the distribution of denudation rates and k sn when using regionally averaged values for model parameters (Fig. 6b). Notably, taking into consideration the variation of runoff and discharge variability for individual sam-ples ( Fig. S5) does not allow identifying significant trends in the data (Adams et al., 2020), as their co-variation with denudation rates mostly reflects their underlying spatial distribution and differences between the margin and the interior of the massif. Even if data for the Cévennes basins appear to be consistent with model predictions for low discharge variability k values, in other areas the k values which would yield the better correspondence between data and the model do not agree with observations. For example, basins from the Ardèche mountains and the Loire headwaters display low and high k values, respectively, whereas the theoretical relationship between denudation rate and k sn in Figs. 6b and S5 would imply the opposite. For some areas from the interior of the massif, we also note high variability in k and runoff for samples presenting similar denudation rates and k sn values (Fig. S5b). Such a discrepancy highlights the limitations of using regionally averaged values for model parameters, such as erodibility or runoff, and notably the necessity of taking into account the strong climatic spatial gradients observed over our study area (Fig. 1b).
Although our dataset does not display a single clear trend between k sn and denudation it is consistent with the global trend delineated by other similar studies using catchmentwide denudation rates (Fig. 9). Our dataset lies at the lower end of the global distribution and overlaps with data from other areas such as the Ecuadorian Andes (Vanacker et al., 2015) and eastern Tibet (Ouimet et al., 2009). We note that the dispersion of our data around the trend, in terms of both denudation rate and k sn , is similar to that observed in other regional studies. As already noted in the compilation of Lague (2014), incision rates present much higher values for a given steepness index range when compared with catchment-wide denudation rates. Such a systematic contrast can arise for various reasons, such as differences between integration timescales of the various types of measurements or  Wobus et al. (2005), Godard et al. (2010Godard et al. ( , 2012Godard et al. ( , 2014, Finnegan et al. (2008), Ouimet et al. (2009), andHarkins et al. (2007), as compiled by Scherler et al. (2017). Additional denudation data from Cyr et al. (2010) are also plotted, as are incision rate data from Avouac (2000, 2001), Kirby and Whipple (2001), and Yanites et al. (2010). a transient landscape response associated with a partial decoupling between channels and hillslopes (e.g., Willenbring et al., 2013;Clubb et al., 2020).

Performance and limitations of modeling formulations
We now discuss the relative behavior of various model formulations of increasing complexity and their ability to explain the variation in denudation rates observed in our dataset. Similarly to results by Campforts et al. (2020), we observe an improvement of the predictive power between the A-SPM and R-SPM types of models ( Fig. 7a and b), highlighting the importance of the spatial distribution of precipitation in modulating the efficiency of fluvial incision processes, at least in this type of context (e.g., Finlayson et al., 2002;Moon et al., 2011;Ferrier et al., 2013a;Adams et al., 2020). Several observations in tectonically active settings (e.g., Godard et al., 2014;Scherler et al., 2014) have advocated for a primary control of rock uplift gradients in dictating the distribution of denudation rates, which is to be expected in a steady-state landscape. Nevertheless, recent results by Adams et al. (2020) also highlight the modulating influence of mean annual precipitation on erodibility and fluvial relief in similar high-uplift settings. Our results support a similar modulation by precipitation in slower-evolution landscapes such as the margin of the Massif Central. In our case, using an ST-SPM setting, the critical shear stress according to observations to τ c = 40 Pa induces a clear failure of the model to reproduce the observed denudation rates (Fig. 7c), and when this parameter is set free the optimization converges toward a very low and actually negligible value (Fig. 7e), which makes this model equivalent to the R-SPM case. This result highlights a clear limitation of the ST-SPM approach when compared with much simpler formulations of the stream power model in the type of relatively low-denudation context we investigate here. It is noteworthy that, in the case of the τ c fixed at 40 Pa, the mismatch appears to be largest for the low-relief catchments of the Loire and Garonne headwaters, while this parameterization still allows predicting denudation rates within ±20 % of the observed uncertainty range for a large proportion of the Cévennes and Ardèche catchments. Field measurements of bedload size show a positive correlation between D 50 and steepness index (Fig. 5b), which would be consistent with smaller incision thresholds in low-relief areas, in contrast with the coarser bedload of the higher-relief Cévennes and Ardèche mountains. Accounting for threshold effects was a critical element in the explanation of the nonlinearity between denudation rates and k sn observed by DiBiase and Whipple (2011), whereas in our case, it appears that a large part of the scatter and deviation from linearity is controlled by the spatial variations in runoff, discharge variability, and erodibility. We note that approximately linear relationship (n ≈ 1) have been reported in other contexts with < 100 mm kyr −1 denudation rates (e.g., Ferrier et al., 2013a;Godard et al., 2019), which might play down the importance of threshold effects in explaining the variability of denudation in such settings. The next increment in model complexity is to account for possible spatial variations in bedload properties. DiBiase et al. (2018) demonstrated a direct impact of contrasts in the size distribution of sediment delivered by hillslopes to channels on landscape relief at two locations in southern California. Attal et al. (2015) documented a power-law relationship between flow competence and D 50 along the Feather River in California, supporting the idea that the intensity of surface processes fluxes, including weathering and erosion, is a primary control on bedload size, and Shobe et al. (2018) explored the theoretical implications of incision thresholds varying with erosion rates. Field measurements in the area we investigate are consistent with such dependency (Fig. 5b), and following Scherler et al. (2017) we have tested the impact of introducing a variable incision threshold dependent on k sn (Eqs. 13 and 3). Such a parameterization yields a significantly better fit than the constant threshold case (Fig. 7d), in particular for some of the low-denudation catchments of the dataset, and allows taking into account differences in bedload size across the divide. However, this model formulation still largely underpredicts erosion rates for many of the low-denudation interior catchments and as a consequence yields a fit to the observations which is still inferior to R-SPM predictions. We note that the stochastic-threshold stream power model is built on the assumption that river incision is detachment-limited and a function of bed shear stress τ . We made many observations of exposed bedrock submitted to active incision processes at our sampling sites, in particular in the high-relief areas east of the divide (Figs. S2 and S4). However, we can consider the possibility that, for some of the low-denudation and lowrelief catchments located west of the divide, a departure from strictly detachment-limited conditions could temporarily occur over the integration timescale of CRN measurements, in particular in relation to post-glacial changes in sediment fluxes. This change in behavior could be an explanation for the failure of the ST-SPM or SVT-SPM to predict observed denudation rates for these catchments unless a very low or negligible incision threshold is considered.
As we have shown earlier, considering k 50 and q (Eq. 13) to be free parameters allows obtaining an adjustment comparable to the R-SPM case. The corresponding D 50 values are lower than present-day observations by a factor of 2, but the ∼ 0.8 exponent of the D 50 − k sn scaling is strikingly similar to observations. It is important to stress that our denudation rates integrate over several tens of thousands of years, a time frame over which important changes in bedload might have occurred, depending on the nature and efficiency of the weathering processes responsible for producing sediment on hillslope (Marshall et al., 2015). Another factor which might contribute to the underprediction of the D 50 is that all the data used here come from surface bedload counts. Therefore, they reflect the surface grain size, which might not be representative of the total volumetric D 50 . Lastly, in order to avoid introducing too many free parameters in the model optimization we have used the same scaling parameters for Eq. (13) over all of our study area. Nevertheless, it could also be noted that the major lithological difference between the Cévennes schists and the granitic-gneissic bedrock of the other regions could have some influence on the spatial variations in bedload size. Indeed, the highly foliated and anisotropic schists produce sediments with a very specific shape factor when compared to the rounded pebbles observed in Ardèche, but such a difference is not resolvable from the available bedload data (Fig. 5b).
As highlighted by Campforts et al. (2020) in the northern Andes, our results confirm the importance of documenting and taking into account spatial variations of environmental parameters. In our case, we observe that gradients in discharge variability and incision threshold, which are key parameters of the stochastic incision family of models, play an important role in explaining the distribution of denudation rates. A key finding of our study is the overall very good performance of R-SPM when compared with more complex formulations of the stream power model.

Importance and definition of threshold effects
Here we discuss the relative importance of incision threshold effects in a broader context by comparing our results with earlier similar studies. In order to carry out this comparison of the different settings in which ST-SPMs have been investigated with similar methodologies as the one we use in our study, we first compute the ratio of erosion rates E to incision threshold (Fig. 10a).
As discussed in Lague et al. (2005) and DiBiase and Whipple (2011) the E/ ratio allows distinguishing three domains defining the importance of threshold effects from domain I with low E/ (where the ST-SPM formulation is relevant and threshold effects are important) to domain III with high E/ (where a simple formulation such as A-SPM or R-SPM adequately describes the incision processes). Domain II is a transitional regime between these two situations. In the case of our study, we distinguish between the two scalings we have used for D 50 calculation ( Fig. 7d and f). We note that E/ varies over 4 orders of magnitude between these different studies. Only our results when using the D 50 scaling set from field observations and the data from DiBiase and Whipple (2011) in the San Gabriel Mountains are clearly in threshold-dominated regime I (Fig. 10a). It appears that all cases in which the incision threshold is a free parameter for model optimization, with either a constant or variable τ c , fall into the transitional regime II. The two cases we consider for Figure 10. (a) Long-term measured erosion rate E normalized by incision threshold (Eq. 7), allowing for distinguishing the various regime domains discussed in Lague et al. (2005) and DiBiase and Whipple (2011). The parameters used in the calculation (k e , τ c ) were obtained from each of the corresponding studies (DiBiase and Whipple, 2011;Scherler et al., 2017;Campforts et al., 2020). For the Ecuadorian Andes data, we use parameters listed in Campforts et al. (2020, (Fig. 7d) and free (Fig. 7f) D 50 scaling (Eq. 13), respectively. (b) Plot of the discharge variability parameter k against normalized critical discharge Q * c calculated from Eq. (11). Black lines correspond to the return time of normalized critical discharge t r (Q * c ) calculated from Eq. (17) (Lague et al., 2005;Lague, 2014). Purple triangles and grey rectangles correspond to DiBiase and Whipple (2011) andCampforts et al. (2020) results, respectively, with the same model parameters as in the previous panel. Circles and squares correspond to our results with observed and free D 50 scalings, respectively. The results from Scherler et al. (2017) were not plotted here, as they use a different definition of the discharge probability distribution.
our Massif Central dataset illustrate the very strong influence of the hypothesis made for bedload size in influencing threshold effects and the associated regime, as there is a factor of 2 in D 50 for values between the two scalings but a corresponding order of magnitude difference in the E/ ratios. Following Lague (2014) we also consider the implications of the relative influence of incision thresholds on the recurrence time t r for incision events (Fig. 10b) as a function of the critical normalized discharge Q * c = Q c /Q (Eq. 11) and the discharge variability parameter k (Eq. 12), where is the incomplete gamma function. We observe that studies in which the threshold is considered to be a free parameter set by model optimization imply very low t r of less than 1 week, which, in the case of our study area in the Massif Central, is not consistent with observations of a limited number of incision events per year. Using the D 50 scaling calibrated with field data yields t r values of a few months for most basins, which is more realistic given the hydrologic regime of our study area. While measurements of bedload size are relatively straightforward to acquire in the field, definition of the appropriate D 50 value to use in ST-SPMs is complicated by the usual high level of dispersion in observations as illustrated by data for the Massif Central (Fig. 5b) and the San Gabriel Mountains (DiBiase and Whipple, 2011, Fig. 6), as well as the high level of variability in space and time of bedload characteristics. Reliable estimates are thus difficult to obtain for τ c , even though the a = 1.5 exponent in Eq. (1) implies that limited changes in D 50 can have important implications for the magnitude of threshold effects. Some earlier studies have avoided this limitation by considering τ c , or the parameters involved in its calculation, to be free with values defined by the model optimization over their dataset. In the case of a spatially variable τ c , our results illustrate a situation in which using such an optimization approach improves the quality of the fit to the erosion rate data (Fig. 7f vs. Fig.7d) but imply a low threshold and unrealistically short return times of incision events with respect to observed present-day conditions. A key feature of our study is that we attempted to constrain every parameter involved in the models and their spatial variations, in particular for the definition of incision thresholds. The importance of such thresholds appears to be strongly contrasting across the geomorphic and climatic gradient of the SE Massif Central.

Conclusions
Our study contributes to the ongoing debate concerning the impact of the distribution of discharge events on river incision and landscape evolution. We provide a case in addition to the handful of studies which have attempted to confront stochastic-threshold incision models with actual observations by focusing on the singular setting of the SE margin of the Massif Central and its pronounced gradient in discharge variability; we try to provide as many field-based controls as possible on all the environmental factors involved in the models, in particular their spatial variations.
Denudation rates are contrasted across our study area and display clear co-variations with precipitation and relief. In contrast, the relationship between denudation and steepness index does not delineate a single trend, suggesting complex interference between the spatial pattern of lithology and climate. A simple version of the stream power model, accounting for the distribution of precipitation, performs better at reproducing observed denudation rates when compared to the more elaborate stochastic-threshold SPM. The performances of such ST-SPMs are critically impacted by the definition and parameterization of the incision threshold. These results raise questions about the reliability of the assumptions usually made regarding bedload properties in such models and the signification of their estimate from field observations in terms of both underlying spatial variations and time integration of drastic changes in surface processes over the Late Pleistocene.
A key endeavor of quantitative geomorphology research is to provide a physically based explanation for the variability of erosion, denudation, or incision rates observed at the Earth's surface, accounting for landscape topographic attributes and environmental forcing of climatic or tectonic origins. The complexity of processes driving river incision has been extensively documented and explored with functional relationships attempting to capture specific mechanisms. The relative success of simpler representations in explaining landscape denudation in our study more generally suggests that the upscaling of such complex mechanistic models toward larger spatial scales requires further investigation.
Data availability. Geoscientific datasets such as denudation rates, morphometric data, and hydrological data are available in the Supplement.