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 fluvial 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 field 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 significant 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 ¹⁰Be cosmogenic radionuclide (CRN) concentration measurements in river sediments, morphometric analysis of river long profiles, and field 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⁻¹ (equivalent to mm/ka in the figures) 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₅₀, 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 field studies and notably the impact of actual constraints on bedload size.

1 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²), 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 short-wavelength 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.

2 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 τ,

$$\begin{array}{}\text{(1)}& I=k_{\mathrm{e}}\left({\mathit{\tau }}^{\mathrm{a}}-{\mathit{\tau }}_{\mathrm{c}}^{\mathrm{a}}\right),\end{array}$$

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

$$\begin{array}{}\text{(2)}& \mathit{\tau }=k_t{\left({\displaystyle \frac{Q}{W}}\right)}^{\mathit{\alpha }}{S}^{\mathit{\beta }},\end{array}$$

with k_t, α, and β as parameters depending on the formulation used to express flow resistance: for example, $\mathit{\alpha }=\mathrm{3}/\mathrm{5}$ and $\mathit{\beta }=\mathrm{7}/\mathrm{10}$ when using Manning's frictional relationship or $\mathit{\alpha }=\mathit{\beta }=\mathrm{2}/\mathrm{3}$ for a Darcy–Weisbach relation (Howard, 1994).

The critical shear stress can be expressed as

$$\begin{array}{}\text{(3)}& {\mathit{\tau }}_{\mathrm{c}}={\mathit{\tau }}_{\mathrm{c}}^{*}g\left({\mathit{\rho }}_{\mathrm{s}}-{\mathit{\rho }}_{\mathrm{w}}\right)D_{\mathrm{50}},\end{array}$$

where D₅₀ is median bedload grain size, ${\mathit{\tau }}_{\mathrm{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

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{W}{W_{\mathrm{b}}}}={\left({\displaystyle \frac{Q}{\stackrel{\mathrm{‾}}{Q}}}\right)}^{{\mathit{\omega }}_{\mathrm{s}}},\end{array}$$

with ω_s an empirical parameter and W_b the bankfull width, scaling with the mean annual discharge $\stackrel{\mathrm{‾}}{Q}$ as

$$\begin{array}{}\text{(5)}& W_{\mathrm{b}}=k_{\mathrm{w}}{\stackrel{\mathrm{‾}}{Q}}^{{\mathit{\omega }}_{\mathrm{b}}},\end{array}$$

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 $\stackrel{\mathrm{‾}}{Q}$ can be calculated from mean runoff $\stackrel{\mathrm{‾}}{R}$ and drainage area A as $\stackrel{\mathrm{‾}}{Q}=\stackrel{\mathrm{‾}}{R}{A}^{\mathrm{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,

$$\begin{array}{}\text{(6)}& I=k_{\mathrm{e}}{k}_t^{\mathrm{a}}{k}_{\mathrm{w}}^{-\mathrm{a}\mathit{\alpha }}{\stackrel{\mathrm{‾}}{R}}^{\mathrm{a}\mathit{\alpha }\left(\mathrm{1}-{\mathit{\omega }}_{\mathrm{b}}\right)}{\left({\displaystyle \frac{Q}{\stackrel{\mathrm{‾}}{Q}}}\right)}^{\mathrm{a}\mathit{\alpha }\left(\mathrm{1}-{\mathit{\omega }}_{\mathrm{s}}\right)}{A}^{\mathrm{a}\mathit{\alpha }\left(\mathrm{1}-{\mathit{\omega }}_{\mathrm{b}}\right)}{S}^{\mathrm{a}\mathit{\beta }}-k_{\mathrm{e}}{\mathit{\tau }}_{\mathrm{c}}^{\mathrm{a}}.\end{array}$$

This equation can be simplified by defining an area exponent $m=a\mathit{\alpha }(\mathrm{1}-{\mathit{\omega }}_{\mathrm{b}})$, a slope exponent n=aβ, an erodibility term $K=k_{\mathrm{e}}{k}_t{k}_{\mathrm{w}}^{-a\mathit{\alpha }}{\stackrel{\mathrm{‾}}{R}}^m$, $\mathit{\gamma }=a\mathit{\alpha }(\mathrm{1}-{\mathit{\omega }}_{\mathrm{s}})$, and a threshold term $\mathrm{\Psi }=k_{\mathrm{e}}{\mathit{\tau }}_{\mathrm{c}}^{\mathrm{a}}$,

$$\begin{array}{}\text{(7)}& I=K{\left({\displaystyle \frac{Q}{\stackrel{\mathrm{‾}}{Q}}}\right)}^{\mathit{\gamma }}{A}^m{S}^n-\mathrm{\Psi }.\end{array}$$

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_{\mathrm{s}}=A^{m/n}S$ (Kirby and Whipple, 2012), Eq. (7) can also be written as

$$\begin{array}{}\text{(8)}& I=K{\left({\displaystyle \frac{Q}{\stackrel{\mathrm{‾}}{Q}}}\right)}^{\mathit{\gamma }}{k}_{\mathrm{s}}^n-\mathrm{\Psi }.\end{array}$$

Assuming that the threshold term Ψ is negligible and that discharge Q is constant and equal to mean discharge $\stackrel{\mathrm{‾}}{Q}$, Eq. (7) simplifies to the classical stream power form, which is widely used in quantitative geomorphology studies:

$$\begin{array}{}\text{(9)}& I=K{A}^m{S}^n.\end{array}$$

Long-term incision E can be calculated by weighting instantaneous incision (Eqs. 7 or 8) with discharge probability pdf(Q) and integrating over Q,

$$\begin{array}{}\text{(10)}& E=\underset{Q_{\mathrm{c}}}{\overset{Q_{\mathrm{m}}}{\int }}I(Q,S,A)\mathrm{pdf}\left(Q\right)\mathrm{d}Q,\end{array}$$

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):

$$\begin{array}{}\text{(11)}& Q_{\mathrm{c}}=\stackrel{\mathrm{‾}}{Q}{\left({\displaystyle \frac{\mathrm{\Psi }}{K}}\right)}^{\frac{\mathrm{1}}{\mathit{\gamma }}}{A}^{-\frac{m}{\mathit{\gamma }}}{S}^{-\frac{n}{\mathit{\gamma }}}=\stackrel{\mathrm{‾}}{Q}{\left({\displaystyle \frac{\mathrm{\Psi }}{K}}\right)}^{\frac{\mathrm{1}}{\mathit{\gamma }}}{k}_{\mathrm{s}}^{-\frac{n}{\mathit{\gamma }}}.\end{array}$$

In incision models taking into account discharge variability, the probability distribution of discharges is often modeled using an inverse gamma law,

$$\begin{array}{}\text{(12)}& \mathrm{pdf}\left(Q^{*}\right)={\displaystyle \frac{k^{k+\mathrm{1}}}{\mathrm{\Gamma }(k+\mathrm{1})}}\exp \left(-{\displaystyle \frac{k}{Q^{*}}}\right)Q^{*-(k+\mathrm{2})},\end{array}$$

where $Q^{*}=Q/\stackrel{\mathrm{‾}}{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.

3 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).

Figure 1 (a) Topographic map of the southeastern margin of the Massif Central (France), 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.

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).

Figure 2 (a) Swath profile for elevation and mean annual precipitation (MAP) across the southeastern margin of the Massif Central (width 60 km). The location of the profile is indicated in Fig. 1b. Solid lines and envelopes are the medians and interquartile ranges, respectively. (b) Values of the discharge variability parameter at gauging stations (green circles). The black solid line is the median along the same swath profile as the panel above. The envelope denotes the interquartile range. The red line corresponds to the interpolated thin plate spline surface (Fig. 1b).

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⁻¹. 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 ¹⁰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.

4 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.

4.1 Discharge variability

4.1.1 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 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 $\mathrm{pdf}\left(\stackrel{\mathrm{‾}}{Q}\right)\propto {\stackrel{\mathrm{‾}}{Q}}^{-\mathit{\alpha }}$. The value of exponent α is an indicator of the importance of high-discharge 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 $\mathit{\alpha }=k+\mathrm{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 long-wavelength 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.

4.1.2 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).

Figure 3 Normalized daily discharge exceedance frequency for selected gauging stations (dark grey circles, see Fig. 1b for location). Dark red lines are fitted inverse gamma distributions (Eq. 12). The corresponding value of the discharge variability parameter k is indicated in each panel. The vertical dashed lines indicate for each station the lower bound of the normalized discharge values used to fit the power law (Clauset et al., 2007).

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 moving 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, 2020; Malcles 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.

4.2 Runoff

4.2.1 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 basin-averaged precipitation and runoff at these instrumented locations in order to be able to compute runoff at the sites where we obtained ¹⁰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².

4.2.2 Results

We observe a positive correlation between basin-averaged MAP and runoff measured at gauging stations (Fig. 4).

Figure 4 Comparison between runoff and mean annual precipitation for 83 selected gauging stations and their corresponding watersheds. Blue and orange symbols correspond to watersheds located on the southeastern margin and inside the Massif Central, respectively. Grey circles indicate basins draining carbonate bedrock. The size of the circles depends on the area of the watersheds. The black dashed line corresponds to $\mathrm{MAP}=\stackrel{\mathrm{‾}}{R}$, and the black solid line is a linear fit excluding carbonate bedrock watersheds.

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 $\stackrel{\mathrm{‾}}{R}$ in order to compute an estimate for runoff $\stackrel{\mathrm{‾}}{R}$ for the catchments where we obtained ¹⁰Be denudation rates (Table S3).

4.3 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² and discharge > 0.1 m³ s⁻¹.

4.3.2 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).

Figure 5 (a) Bankfull river width against discharge calculated using the calibration in 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₅₀ 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.

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=\mathrm{0.45}$. In this case we obtain a value of $k_{\mathrm{w}}=\mathrm{5.56}\pm \mathrm{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.

4.4 River long profile

4.4.1 Methods

We extracted the catchment contour and upstream river network for each site where ¹⁰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² 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 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).

4.4.2 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.

4.5 Bedload characteristics

4.5.1 Methods

We use bedload median size D₅₀ 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₅₀ 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₅₀ 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₅₀ (Attal et al., 2015; Scherler et al., 2017), such as

$$\begin{array}{}\text{(13)}& D_{\mathrm{50}}=k_{\mathrm{50}}{k}_{\mathrm{sn}}^q.\end{array}$$

4.5.2 Results

Median bedload size ranges from 10 to 200 mm over our study area. The average D₅₀ 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₅₀ and k_sn, with exponent ∼0.8. We do not observe a difference in scaling between the two sides of the divide. For the low-k_sn regions, corresponding mostly to areas west of the divide, measured D₅₀ values are dominantly in the 30–60 mm range, whereas they are mostly >50 mm for the eastern catchments draining toward the Rhône.

Figure 6 (a) Comparison between erosion rate and basin slope. Symbols are colored according to the location of the measurements (Cévennes, Ardèche, Loire, and Garonne). Basin CDX-30 is a very small catchment when compared to the rest of our dataset and is colored in yellow. (b) Comparison between normalized channel steepness index and denudation rates. Black lines correspond to the theoretical relationship between denudation rate and steepness index computed from Eqs. (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_{\mathrm{e}}=\mathrm{8}\times {\mathrm{10}}^{-\mathrm{13}}$ m^2.5 s² kg^−1.5. We use a regionally averaged value for runoff of $\stackrel{\mathrm{‾}}{R}=\mathrm{650}$ mm a⁻¹.

4.6 Cosmogenic nuclides

4.6.1 Methods

We measured ¹⁰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 basins were primarily selected according to their size (21 to 95 km² 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 ⁹Be carrier solution, and Be was extracted with ions exchange columns. ${}^{\mathrm{10}}\mathrm{Be}{/}^{\mathrm{9}}\mathrm{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.

4.6.2 Results

Measured ¹⁰Be concentrations in the 34 sampled catchments range from 51 to 345×10³ at g⁻¹ (Table S1 in the Supplement). The corresponding calculated denudation rates range from 24 to 126 mm kyr⁻¹ (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⁻¹, respectively. Significantly faster denudation is observed along the margin of the massif, in the Ardèche mountains (58±7 mm kyr⁻¹), and especially in the Cévennes area (95±11 mm kyr⁻¹). 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, 2002; Molliex et al., 2016; Olivetti et al., 2016).

Denudation rates are positively correlated with basin-average 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², 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).

Table 1 Reference parameters used in the modeling.

Note that ${\mathit{\theta }}_{\mathrm{ref}}\equiv m/n=\mathit{\alpha }\left(\mathrm{1}-{\mathit{\omega }}_{\mathrm{b}}\right)/\mathit{\beta }$.

5 Denudation rate modeling

5.1 Modeling approach

We model denudation rates for each basin where we have ¹⁰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

$$\begin{array}{}\text{(14)}& E=K{k}_{\mathrm{s}}^{\mathrm{n}}=k_{\mathrm{e}}{k}_t{k}_{\mathrm{w}}^{-a\mathit{\alpha }}{\stackrel{\mathrm{‾}}{R}}^m{k}_{\mathrm{s}}^{\mathrm{n}}.\end{array}$$

In this case we can either neglect the spatial variations of runoff, consider $\stackrel{\mathrm{‾}}{R}$ to be constant, and set it to a regional average 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₅₀ and its relationship with k_sn, and we first set the parameters k₅₀ and q of Eq. (13) according to these constraints (Fig. 5b). In a second run the models consider k₅₀ 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² for a linear fit between O and M or the χ² value,

$$\begin{array}{}\text{(15)}& {\mathit{\chi }}^{\mathrm{2}}=\sum _{i=\mathrm{1}}^n{\left({\displaystyle \frac{O_i-M_i}{{\mathit{\sigma }}_i}}\right)}^{\mathrm{2}},\end{array}$$

and we maximize Nash–Sutcliffe (NS) model efficiency for our determination of optimal parameters,

$$\begin{array}{}\text{(16)}& \mathrm{NS}=\mathrm{1}-{\displaystyle \frac{\sum _{i=\mathrm{1}}^n{\left(O_i-M_i\right)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^n{\left(O_i-\stackrel{\mathrm{‾}}{O}\right)}^{\mathrm{2}}}}.\end{array}$$

5.2 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⁻¹).

Figure 7 Comparison between measured and modeled denudation rates for various formulations of the stream power model. Black lines correspond to a perfect match between observed and modeled denudation rates, and dashed lines indicate ±20 % and ±50 % deviations. Basin 30, with a much smaller area than the rest of the dataset, is outlined in yellow. (a) Area-based stream power model (A-SPM). The erodibility coefficient k_e and modulation factor for Cévennes schists f_s are free parameters. (b) Runoff-based stream power model (R-SPM). (c) Stochastic stream power model with an imposed spatially constant threshold τ_c (ST-SPM) defined using the mean regional average D₅₀ value (Eq. 3 and Fig. 5b). (d) Stochastic stream power model, with a spatially variable threshold (SVT-SPM), scaled with k_sn values (Eqs. 13 and 3). The scaling parameters are set according to the observed relationship between D₅₀ and k_sn values in Fig. 5b. (e) Stochastic stream power model (same as c), but with τ_c as a free parameter (in addition to k_e and f_s). (f) Stochastic stream power model, with a spatially variable threshold (SVT-SPM, same as d) but with the scaling parameters k₅₀ and q from Eq. (13) being free (in addition to k_e and f_s).

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₅₀ (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₅₀ and k_sn of ∼0.8 is almost identical to the observed value (Fig. 5b), while the k₅₀ factor implies ∼2 times lower predicted D₅₀ values, but which are still within the lower range of observations.

6 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).

6.1 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⁻¹ 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⁻¹, 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⁻¹ 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²) 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.

Figure 8 (a) Denudation rates against basin mean annual runoff calculated using the calibration of Fig. 4. Symbols are colored according to the four different subregions, and the black solid line is a linear fit to the data. (b) Basin relief against mean annual runoff. Symbols are colored according to the denudation rates of the basins.

Denudation rates are positively correlated with basin-average 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 samples (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).

Figure 9 Comparison between our results and other published datasets for normalized channel steepness index (θ_ref=0.45) and basin-wide denudation rate calculated from ¹⁰Be concentrations. Incision rates calculated from terraces in different climatic and tectonic contexts are also plotted (green triangles). The black solid line is a power-law fit to the aggregated basin-wide denudation rate dataset. Original data for the Ecuadorian Andes are from Vanacker et al. (2015) as reported in Campforts et al. (2020). Original data for the San Gabriel Mountains are from DiBiase et al. (2010) as reported in DiBiase and Whipple (2011). Original data for the Himalayas and eastern Tibet are from Scherler et al. (2014), Wobus et al. (2005), Godard et al. (2010, 2012, 2014), Finnegan et al. (2008), Ouimet et al. (2009), and Harkins 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 Lavé and Avouac (2000, 2001), Kirby and Whipple (2001), and Yanites et al. (2010).

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 catchment-wide 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 a transient landscape response associated with a partial decoupling between channels and hillslopes (e.g., Willenbring et al., 2013; Clubb et al., 2020).

6.2 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₅₀ 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⁻¹ 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₅₀ 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 low-relief 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₅₀ and q (Eq. 13) to be free parameters allows obtaining an adjustment comparable to the R-SPM case. The corresponding D₅₀ values are lower than present-day observations by a factor of 2, but the ∼0.8 exponent of the D₅₀−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₅₀ 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₅₀. 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.

6.3 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).

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, Table 4) for their variable runoff, variability, and erodibility ST-SPM model (model 8). For the Himalayas and eastern Tibet data, we use the parameters listed in Scherler et al. (2017, Table 5) for their “Himalaya and eastern Tibet TRMM models”, with either constant D₅₀ and constant ${\mathit{\tau }}_{\mathrm{c}}^{*}$ (solid black line) or variable D₅₀ and constant ${\mathit{\tau }}_{\mathrm{c}}^{*}$ (dashed black line). Dark and light green curves correspond to our sampled basins with observed (Fig. 7d) and free (Fig. 7f) D₅₀ scaling (Eq. 13), respectively. (b) Plot of the discharge variability parameter k against normalized critical discharge $Q_{\mathrm{c}}^{*}$ calculated from Eq. (11). Black lines correspond to the return time of normalized critical discharge $t_{\mathrm{r}}\left(Q_{\mathrm{c}}^{*}\right)$ calculated from Eq. (17) (Lague et al., 2005; Lague, 2014). Purple triangles and grey rectangles correspond to DiBiase and Whipple (2011) and Campforts 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₅₀ scalings, respectively. The results from Scherler et al. (2017) were not plotted here, as they use a different definition of the discharge probability distribution.

As discussed in Lague et al. (2005) and DiBiase and Whipple (2011) the $E/\mathrm{\Psi }$ ratio allows distinguishing three domains defining the importance of threshold effects from domain I with low $E/\mathrm{\Psi }$ (where the ST-SPM formulation is relevant and threshold effects are important) to domain III with high $E/\mathrm{\Psi }$ (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₅₀ calculation (Fig. 7d and f). We note that $E/\mathrm{\Psi }$ varies over 4 orders of magnitude between these different studies. Only our results when using the D₅₀ 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 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₅₀ for values between the two scalings but a corresponding order of magnitude difference in the $E/\mathrm{\Psi }$ 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_{\mathrm{c}}^{*}=Q_{\mathrm{c}}/\stackrel{\mathrm{‾}}{Q}$ (Eq. 11) and the discharge variability parameter k (Eq. 12),

$$\begin{array}{}\text{(17)}& t_{\mathrm{r}}\left(Q_{\mathrm{c}}^{*}\right)=\mathrm{\Gamma }{\left(k/Q_{\mathrm{c}}^{*},k+\mathrm{1}\right)}^{-\mathrm{1}},\end{array}$$

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₅₀ 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₅₀ 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₅₀ 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.

7 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.