CC-BY 4.0
1. Introduction
Landscapes evolve primarily through the interplay of tectonics and climatic forcings, which drive surface processes such as weathering, erosion, and sediment transport (Champagnac et al., 2012; Kirby and Whipple, 2012). While many theoretical formulations conceptualize surface processes using a steady-state framework, where these forcings balance each other, widespread evidence points to ubiquitous transient changes in surface-process intensity and topographic properties across natural landscapes (Mudd, 2016; Prince and Spotila, 2013; Godard, Salgado, et al., 2021). This unsteadiness manifests across all spatial scales, from the response of individual hillslopes to the reorganization of entire drainage systems, and across a range of timescales spanning hundreds to millions of years. Understanding these transient responses is critical in particular over Late Quaternary timescales, during which recent climatic changes—notably associated with the Pleistocene–Holocene transition—have strongly influenced landscape dynamics and set the stage for the geomorphological systems as we observe them today (Mariotti et al., 2021; Madoff and Putkonen, 2022). Nevertheless, the precise nature and strength of the connection between climate change and erosion remain subjects of active debate, hindering the analysis of recent landscape evolution and its interactions with ecosystems and societies.
To characterize these transient erosion responses and understand landscape dynamics across timescales, researchers often rely on the analysis of sedimentary archives. These records form through the progressive accumulation of erosion products and carry a wealth of information on the conditions and processes acting in the contributing catchment at various times (e.g. Kober et al., 2019; Mariotti et al., 2021). However, several limitations exist: incomplete or disturbed records, sparse sampling density, poor control on sediment provenance, and the sometimes significant distance between sediment source area and depositional sink. Furthermore, interpreting these records can be challenging when deconvolving the distinct contributions of tectonic, climatic and anthropogenic influences, which are often mixed and blurred. These constraints underscore the need to develop complementary approaches, especially those capable of documenting near-field erosion dynamics, as close as possible to the active sites of surface processes (Godard, Mudd, et al., 2025).
Terrestrial Cosmogenic Nuclides (TCN) have emerged as a highly effective tool for investigating the dynamics of surface processes, enabling quantitative estimations of past and present denudation rates (von Blanckenburg, 2005; Granger and Schaller, 2014; Schaefer et al., 2022; Codilean et al., 2022). While a single TCN measurement only provides a single denudation rate, the use of multiple TCN with distinct radioactive decay rates offers invaluable insights into temporal changes in surface processes (e.g. Hippe, Kober, et al., 2012; Schimmelpfennig et al., 2022; Godard, Siame, et al., 2024). Among these, paired measurements of the terrestrial cosmogenic nuclides 10Be and 14Cin situ, both produced in situ within the quartz mineral, are particularly promising for resolving Late Pleistocene and Holocene landscape dynamics, over 1 to 10 ka timescales (Mudd, 2016; Hippe, 2017; Skov et al., 2019). The very distinct half-lives of cosmogenic nuclides like 14C (5.7 ka) and 10Be (1.39 Ma) make them ideally suited for quantifying changes in denudation rates across these critical timescales (Fülöp, Bishop, et al., 2015; Hippe, Jansen, et al., 2021). Indeed, the much longer-lived 10Be is sensitive to landscape dynamics on the 10–100 ka timescale, whereas the comparatively short-lived 14Cin situ is mostly sensitive to recent events on the 100 a–10 ka timescale, with no memory of older denudation conditions. This pair is in particular sensitive to recent perturbations, such as land cover changes, landslides, or fluctuations in glacial extent, which can significantly shift the relative concentrations of the two nuclides (14Cin situ/10Be ratio) away from its steady-state equilibrium value (Hippe, 2017; Schimmelpfennig et al., 2022). Despite its considerable potential, very few datasets using the 10Be–14Cin situ system to constrain changes in denudation dynamics have been published to date, primarily because of the limited number of specialized 14Cin situ extraction facilities (Lupker et al., 2019; Fülöp, Fink, et al., 2019; Lifton et al., 2023), as well as the more complex interpretations required by this system (Skov et al., 2019; Hippe, 2017).
The Pleistocene–Holocene transition represents one of the most dramatic recent environmental shifts, profoundly shaping the Earth’s contemporary landscapes and ecosystems (e.g. Shakun and Carlson, 2010; Malhi et al., 2016; Ellis, 2021). During this period, substantial increases in temperature and modifications in precipitation regimes led to significant changes in sediment production and transport mechanisms (e.g. Hughes et al., 2009; Godard, Mudd, et al., 2025). Divergences in the apparent denudation rates derived from 10Be and 14Cin situ can serve as a diagnostic signal for erosional transience during the Holocene, including periods of significant anthropogenic influence. Over the Late Pleistocene and Holocene, landscapes experienced a conjunction of environmental and cultural transformations, with climatic variability often coinciding with the emergence, expansion and collapse of agropastoral societies (e.g. Ruddiman, 2007; Kirch, 2005). In many early agricultural centers, accelerated soil thinning and erosion have been directly linked to human activities such as widespread burning, deforestation, intensive grazing, or cultivation, often exacerbated by changing climatic conditions (e.g. Montgomery, 2007; Owens, 2020). This entanglement of factors poses a significant challenge to accurately discerning the relative contributions of natural environmental change versus anthropogenic forcing to long-term soil sustainability and landscape evolution (Hippe, Jansen, et al., 2021). The idea of an “Early Anthropocene” emphasizes the human and climate-driven impacts on pre-industrial landscapes and highlights this complex and often inseparable relationship (ibid.).
To effectively address these questions, a careful selection of study areas is required. Regions characterized by slow tectonic deformation and low background denudation rates offer optimal conditions for isolating, on 1–10 ka timescales, the influence of climatic variations on surface processes from often overwhelming tectonic signals. In this context, the southeastern French Massif Central (FMC) presents an interesting natural laboratory (Olivetti, Godard, et al., 2016; Desormeaux et al., 2022). This mountainous region features a strong geomorphological gradient from its incised eastern margin to the more subdued relief of the Massif interior. Focusing on small catchments in this type of context allows us to minimize the complexities associated with extensive sediment storage and long transport lags that often complicate the interpretation of detrital TCN signals.
In this study, we present a new set of 14Cin situ measurements in river sediment combined with existing 10Be concentrations from 10 catchments on the eastern margin of the French Massif Central. We analyze the denudation rates derived from these nuclide concentrations and evaluate them against different surface-process scenarios to explain the variations in disequilibrium observed across the margin. We explore in particular the impact of instantaneous variations of denudation rates and soil thickness.
2. Settings
The French Massif Central (FMC) belongs to a series of low- to middle-elevation mountain ranges that represent the remnants of the Hercynian orogeny across western and central Europe, such as the Vosges or Bohemian Massif (Olivetti, Balestrieri, et al., 2020). Due to significant uplift and volcanic activity during the Late Cenozoic (Séranne et al., 2002; Malcles et al., 2020), the FMC stands as a major topographic feature, and is characterized by a very steep southeastern margin with topography rising from a few hundred metres a.s.l. in the Rhône river valley, to more than 1000 m asl in the interior of the FMC (Figure 1).
Location map of the studied basins along the SE margin of the French Massif Central. Basins are colored according to their average slope. The three classes also correspond to the three clusters that were defined independently from cosmogenic-nuclide data, as explained in the Methods. The numbers used for identification are the same as in Desormeaux et al. (2022), where the 10Be data were published. The red line indicates the drainage divide between the Rhône basin to the east, and the Loire and Garonne basins to the west.
The topographic gradient is associated with a morphological contrast between the Ardèche and Cévennes Mountains with steep near-threshold hillslopes and active river incision into bedrock, and the gently rolling hillslopes observed farther west across the margin edge, in the Margeride Mountains (Desormeaux et al., 2022). This steep margin acts as a major topographic barrier for moisture transport from the Mediterranean, and an orographic effect leads to contrasting mean annual precipitation: >1500 mm at the margin and ∼500 mm in the interior of the Massif. The area is also subject to intense storm events several times per year, leading to important discharge variability along the streams that drain the margin (Le Lay and Saulnier, 2007; Desormeaux et al., 2022). Denudation rates derived from 10Be concentrations in river sediments also highlight a strong spatial gradient across the margin. The Ardèche and Cévennes mountains are characterized by denudation rates in the 80–120 mm/ka range, whereas moving west beyond the Rhone catchment divide, denudation rates are usually <40 mm/ka (Olivetti, Godard, et al., 2016; Desormeaux et al., 2022).
Like other middle-elevation Hercynian mountain domains in Europe, the French Massif Central has been subject to limited glaciation during the Late Quaternary, although several groups of glaciers and icefields have been identified, notably in the central part of the FMC, in the Cantal, Mont Dore and Aubrac mountains (Ancrenaz et al., 2022). Further east, closer to the studied area, smaller glacier systems were also recognized in the Margeride and Lozère Mountains, with a limited geomorphic imprint, both in terms of cumulative erosion and associated deposits. This is one of the key advantages of the study area, as the TCN signal is unlikely to be strongly affected by glacial shielding, glacial erosion or massive paraglacial sediment reworking in the investigated basins. This limited influence of the Quaternary Glaciations allows for a focused investigation into the impacts of climatic and anthropogenic factors on weathering and hillslope processes.
3. Methods
3.1. Analytical approach
Desormeaux et al. (2022) sampled 34 medium-sized basins along the southeastern margin of the FMC and measured 10Be concentrations. We selected 10 of these basins for 14Cin situ analysis to investigate their recent dynamics. The basins are affected by a large spectrum of surface processes, due to varying mountainous landforms ranging from gently rolling hillslopes of the FMC interior to the steep margin of the Ardèche and Cévennes Mountains. During the sample selection, Desormeaux et al. (ibid.) took care to avoid areas impacted by active landslides or dense human activity.
In situ 14C extractions were performed at ETH Zürich on 5–6 g aliquots of the 10 quartz samples (250–1000 μm fraction) used for the 10Be analyses (ibid.). The extraction system and method outlined in Lupker et al. (2019) allows degassing of all carbon species, including cosmogenic 14C, from the quartz at temperatures up to 1670 °C. After purification and quantification of the CO2, the 14C/12C ratios were measured at the gas ion source of the MICADAS accelerator mass spectrometer at ETH Zürich (Synal et al., 2007; Fahrni et al., 2013). The 14C concentrations were calculated following the method of Hippe and Lifton (2014) (Table 1). Measurement results of the inter-comparison quartz materials CRONUS-A and CRONUS-N from the ETH extraction line used in this study are reported in Lupker et al. (2019).
14Cin situ analytical results (see text for details)
| Sample | Quartz mass (g) | C Yield (μgC/g) | 𝛿13C (‰) | Fraction modern | 14C (×106 at.) | 14C - bk corr. ( × 106 at.) | [14C] ( × 103 at/g) |
|---|---|---|---|---|---|---|---|
| CDX-CRN-02 | 5.7463 | 2.2 | −5.9 | 0.3035 ± 0.0030 | 0.5400 ± 0.0062 | 0.5240 ± 0.0078 | 91.181 ± 1.353 |
| CDX-CRN-04 | 6.9030 | 13.8 | −13.9 | 0.0386 ± 0.0010 | 0.2545 ± 0.0091 | 0.2384 ± 0.0103 | 34.542 ± 1.487 |
| CDX-CRN-05 | 6.0447 | 20.9 | −9.9 | 0.0357 ± 0.0010 | 0.3033 ± 0.0108 | 0.2872 ± 0.0118 | 47.516 ± 1.951 |
| CDX-CRN-06 | 5.9356 | 58.7 | −20.3 | 0.0158 ± 0.0010 | 0.3352 ± 0.0167 | 0.3191 ± 0.0173 | 53.762 ± 2.921 |
| CDX-CRN-08 | 4.5762 | 1.1 | −5.7 | 0.4728 ± 0.0000 | 0.6195 ± 0.0000 | 0.6034 ± 0.0047 | 131.859 ± 1.023 |
| CDX-CRN-09 | 5.9672 | 3.6 | −9.9 | 0.6048 ± 0.0050 | 1.3857 ± 0.0118 | 1.3697 ± 0.0127 | 229.532 ± 2.124 |
| CDX-CRN-12 | 6.5484 | 3.4 | −7.7 | 0.5876 ± 0.0050 | 1.3795 ± 0.0120 | 1.3635 ± 0.0129 | 208.213 ± 1.967 |
| CDX-CRN-15 | 4.3083 | 94.8 | −22.7 | 0.0065 ± 0.0000 | 0.1598 ± 0.0000 | 0.1437 ± 0.0047 | 33.359 ± 1.087 |
| CDX-CRN-17 | 6.2913 | 3.6 | −9.3 | 0.3025 ± 0.0040 | 0.7188 ± 0.0084 | 0.7028 ± 0.0096 | 111.705 ± 1.530 |
| CDX-CRN-19 | 4.6725 | 3.0 | −10.3 | 0.1600 ± 0.0000 | 0.2959 ± 0.0000 | 0.2799 ± 0.0047 | 59.893 ± 1.002 |
Note that on all figures the basins are labeled with the number at the end of the sample name reported here. All indicated uncertainties are ±1𝜎. The process blank associated with this series of sample contained 16038 ± 4682 atoms of 14C. The subtracted blank was averaged over 20 blank measurements obtained between 06.01.2020 and 27.04.2021. 10Be results and basin properties are reported in Desormeaux et al. (2022).
3.2. Steady-state interpretations of concentrations
Steady-state denudation corresponds to the situation in which a constant surface denudation rate has been applied long enough that the TCN concentration depth profile no longer varies with time. The measured TCN concentrations C (atoms/M), for 10Be and 14Cin situ, can be interpreted in terms of steady-state denudation rates 𝜀 (L/T),
| \begin {equation}\label {eq:ss} C=\sum _i \frac {P_i}{\frac {\rho \varepsilon }{\Lambda _i}+\lambda }, \end {equation} | (1) |
Steady-state apparent denudation rates were calculated using the Lal-Stone time-independent scaling scheme (Stone, 2000), with the calibration dataset of Borchers et al. (2016) for spallation production rates and muon production properties from Braucher, Bourlès, et al. (2013), Braucher, Brown, et al. (2003). We use 160 g/cm2 for the effective neutron attenuation length in rock and a rock density of 2.65 g/cm3. No shielding correction was considered.
To investigate the influence of soil mixing dynamics, we also use the formalism proposed by Foster et al. (2015) for a steady-state mixed soil of thickness H (L), where concentration C (atoms/M) can be calculated as,
| \begin {equation}\label {eq:foster} C = \frac {\frac {\bar {P}H}{\varepsilon \beta } + C_b }{1+\frac {\lambda H}{\varepsilon \beta }}. \end {equation} | (2) |
3.3. Hilltop curvature calculation
Morphological data for the studied basins are reported by Desormeaux et al. (2022). In addition, we computed basin-averaged hilltop curvature CHT (1/L), using open-access 2-m resolution lidar Digital Elevation Models (LIDAR HD; IGN, 2026). Obtaining basin-averaged estimates of hilltop curvature allows us to use the classical relationship in transport-limited hillslopes between CHT and denudation rate 𝜀 (Hurst et al., 2012).
| \begin {equation}\label {eq:cht} \varepsilon = \frac {K |C_{\mathrm {HT}}|}{\beta }. \end {equation} | (3) |
This relationship involves the hillslope transport coefficient K (L2/T), for which published compilations exist (Richardson et al., 2019), and will be used to analyze our dataset.
3.4. Transient dynamics modelling
To study the possible disequilibrium associated with deviations from steady state in the 10Be–14Cin situ system, we also analyze the transient evolution of the concentrations in response to perturbations and how they evolve from one steady-state value toward another. We perform an inversion of the 10Be and 14Cin situ concentrations using a Markov chain Monte Carlo (MCMC) sampler (Vrugt et al., 2009; Hartig et al., 2023). For each model we ran 8 chains in parallel, each 5 × 105 steps in length, with a burn-in phase corresponding to the first 10% of each chain. Convergence was checked using the multivariate potential scale reduction factor (Gelman and Rubin, 1992). The likelihood $\mathcal{L}$ was defined as,
| \begin {equation} \ln \mathcal {L} = \sum _{i} \ln \left ( \frac {1}{\sqrt {2\uppi } \, \sigma _i} \right ) - \frac {1}{2} \sum _{i} \frac {\left (C_i^{\mathrm {obs}} - C_i^{\mathrm {mod}}\right )^2}{\sigma _i^2}, \end {equation} | (4) |
The inversions are constructed according to two distinct transient scenarios inspired by those investigated in Hippe, Jansen, et al. (2021). The first scenario considers a step change in denudation rate at a fixed time T, synchronous for all basins. Contrary to Hippe, Jansen, et al. (ibid.), we formulate the evolution in terms of changes of the hillslope transport coefficient K, and calculate denudation rates using Equation (3). The hillslope transport coefficient is setting the pace of hillslope morphological evolution and is modulated by climatic conditions and vegetation (Richardson et al., 2019). Following Hippe, Jansen, et al. (2021), we consider a two-step scenario. In contrast to their approach, in which the change is specified directly in terms of denudation rate, we instead impose a change in K from K1 to K2 at time T.
The coefficient K is considered to be the same for all basins within each of the clusters, which we define with k-means clustering according to the values of denudation rates derived from the 10Be and 14Cin situ concentrations. We assume that basins in each cluster share similar lithologies and climatic conditions. This approach has the advantage of reducing the dimension of the problem to be solved. For example, in the case of a pair of basins, for which we have determined the 10Be and 14Cin situ concentrations, we have a total of four concentrations, in other words two concentrations per basin. Formulating the problem as a step change in denudation rates at time T requires solving for five unknowns: the transition time T, the denudation rate before the transition for each basin, and the denudation rate after the transition for each basin. The problem is under-determined in this case. Conversely, formulating the evolution in terms of a step change in K reduces the problem to three unknowns: T, K1 and K2, because K is a regional, climatically-controlled parameter whose value at a given time is the same for both basins. We use an exponential prior on the K parameters, based on the compilation of Richardson et al. (2019). The forward model for the prediction of concentrations accounts for time variations in production rates following the Lal-Stone time-dependent scaling scheme (Lal, 1991; Stone, 2000; Balco et al., 2008), and the virtual dipole moment reconstruction of Muscheler et al. (2005).
The second scenario considers an instantaneous soil removal event at a time T, again synchronously for all basins. The thickness of soil removed is identical for basins belonging to the same cluster, as is the hillslope transport coefficient K, which does not change through time.
4. Results
4.1. Analysis of concentrations and steady-state predictions
The measured 14Cin situ concentrations range from 33 ± 1 × 103 to 230 ± 2 × 103 at/g of quartz (Table 1). The lowest concentration is at least one order of magnitude above the +1𝜎 analytical blank value. When combined with the 10Be data of Desormeaux et al. (2022), all our samples plot below the steady-state denudation curve on a two-nuclide plot (solid line, Figure 2A). This curve corresponds to the predicted surface concentrations in the case of a constant denudation through time. This position indicates that the measured TCN concentrations cannot be explained by a steady-state denudation history nor by a constant exposure of samples at the surface. We observe that the amplitude of deviation from the steady-state denudation curve is highly variable, and samples do not cluster as a single group. Some samples are close to steady-state denudation conditions, whereas most samples are distant from the steady-state curve and appear to be very depleted in 14Cin situ. Most of these samples with very low 14Cin situ/10Be are located east of the drainage divide, where topographic gradients are the highest (Figure 1). We computed steady-state denudation rates based on the measured 14Cin situ and 10Be concentrations (Figure 2B). These steady-state rates for the two nuclides are clearly different in most cases, with up to >5-fold differences, the 14Cin situ-derived denudation rates being systematically higher than their 10Be equivalent. This difference can be considered as a clue for an intensification of surface processes in the recent past, within the sensitivity time window of 14C.
(A) Normalized two-nuclide plot for the 10Be–14Cin situ system. The concentrations are normalized according to the individual production rates of the basins. Samples are colored according to the denudation-derived cluster they belong to. Ellipses represent 2𝜎 confidence intervals. The solid and dashed curves correspond to the prediction for steady-state denudation or constant exposure of the surface, that is absence of denudation. (B) Comparison of steady-state denudation rates derived from the 10Be and 14Cin situ measured concentrations (error bars are ±2𝜎). The solid line indicates equal 10Be and 14Cin situ denudation rates, and dashed lines correspond to 14Cin situ denudation rates higher by factors of 2, 5 and 10.
The k-mean clustering based on the apparent denudation rates (Figure 2B) yielded 3 clusters. Cluster 1 corresponds to basins located in the interior of the FMC with low 10Be-derived denudation rates and a moderate disequilibrium between 10Be and 14Cin situ. Cluster 2 corresponds to basins on the topographic margin, mostly in the Cévennes east of the divide, with both high 10Be-derived denudation rates and the strongest disequilibrium. Cluster 3 corresponds to basins on the topographic margin and east of the divide with high 10Be-derived denudation rates and a lower 10Be–14Cin situ disequilibrium. Even though the average slope of the basin was not used for the TCN-based clustering, the three resulting clusters are consistent with the slope classes used in Figure 1 and correspond to a geomorphic gradient of increasing slopes.
We observe a clear correlation (R2 = 0.82) between 10Be-derived denudation rates and hilltop curvature (CHT), except for basin 6, which has the highest curvature (Figure 3). This correlation suggests the importance of transport-limited hillslope dynamics over long time scales, and the obtained transport coefficient K = 0.009 ± 0.001 m2/a is within the range reported by Richardson et al. (2019) for comparable settings. Basin 6 is also the steepest basin, and its deviation from the observed linear trend may result from a transition from transport- to weathering-limited conditions on its hillslopes (Godard, Ollivier, et al., 2016).
Relationship between 10Be steady-state denudation rate (error bars are ±2𝜎) and the absolute value of hilltop curvature extracted from a LiDAR DEM, for the studied basins. Black line is a linear regression (R2 = 0.82 and p < 0.001) and the envelope denotes its 95% confidence interval. Basin 6 was excluded from the regression (justification in text). The slope of the regression corresponds to a hillslope transport coefficient of 0.009 ± 0.001 m2/a.
As a single, basin-wide steady-state denudation rate over the period of resolution cannot explain the concentrations observed for each of our paired 10Be–14Cin situ samples, we explore other processes. Specifically, we investigate steady-state soil mixing and associated dynamic burial, which has been shown to induce apparent disequilibrium in other paired nuclide systems (Makhubela et al., 2019; Knudsen et al., 2019; Hippe, Jansen, et al., 2021; Godard, Siame, et al., 2024). For each basin, we use the formalism proposed by Foster et al. (2015) to compute the denudation rates and soil-mixing thicknesses that explain both the observed 10Be and 14Cin situ concentrations (Figure 4). The modeled denudation rates range from 20 to 120 mm/ka, with mixing thicknesses reaching more than 10 meters for some basins. We did not observe such thick soils in the studied area, even in the low relief interior of the FMC. Furthermore, the fact that the thickest predicted soils are associated with the steepest basins (Figure 4B) suggests that this process of soil mixing is unlikely to be the main factor explaining the observed disequilibrium.
(A) Relationship between steady-state soil mixing depth and denudation rate computed according to the formalism of Foster et al. (2015), using both 10Be and 14Cin situ data. Ellipses correspond to 2𝜎 confidence intervals. (B) Relationship between the steady-state mixing depth and basin average slope (error bars are ±2𝜎).
We also use a simple mixing budget to test the idea that the 14Cin situ–10Be disequilibrium results from a dilution by the tapping of sources of old (>50 ka) and shielded sediments (Figure 5). Deposits such as deep fluvial, glacial or lacustrine fills, can have very low 14Cin situ concentrations and their mixing with the active sediments transiting in the river network could substantially affect the measured concentration. For basins where the disequilibrium is low, a moderate amount (on the order of 10%) of material derived from old and shielded sediments could explain the observed 14Cin situ/10Be ratios. However, for many catchments, the estimated amount is higher than 50%, and no field evidence indicates the existence of reservoirs large enough to be a source of such sediments. We also note that the proportion of these old sediments is highest in the steepest catchments, where intermediate storage of material is least likely. Finally, the existence of such large reservoirs of completely shielded sediments is also highly unlikely from a geomorphological point of view, as any transient deposit would immediately start to accumulate cosmogenic nuclides over depths of several meters.
(A) Steady-state denudation rate ratios between 10Be and 14Cin situ as a function of the estimated fraction of total flux from old sediments. (B) Fraction of the total flux from old sediments as a function of the basin slope.
4.2. Transient dynamics modelling
As steady-state solutions do not provide acceptable explanations for the observed concentrations, we investigate the impact of transient perturbations. As suggested by previous studies (Hippe, Jansen, et al., 2021), the low 14Cin situ/10Be ratios could be explained by a transient response to recent changes in surface processes, either through an increase in denudation rates or through a discrete soil-removal event (ibid.). Indeed, both hypotheses contribute to rapidly bringing rocks with lower 14Cin situ/10Be ratios to the surface, and this disequilibrium can persist for hundreds to thousands of years (Figure 6).
Examples of two possible transient scenarios explaining the observed 10Be and 14Cin situ concentrations for two samples in our dataset (CDX-CRN-02 and CDX-CRN-19). Both scenarios correspond to a perturbation in the past, starting from steady-state conditions. Scenario 1 (orange curves) considers an instantaneous increase in denudation rate, driven by an increase of the hillslope transport coefficient K. Scenario 2 (pink curves) considers an instantaneous removal of a given thickness of soil with no change in the hillslope transport coefficient. The parameters used here are for illustrative purpose and do not correspond to best fitting values from an optimization procedure. The ellipse color corresponds to the cluster the sample belongs to.
In Scenario 1, we consider the implications of a step change in denudation rate at a given transition time T, occurring synchronously in all basins. As explained above, we do not express this change directly in terms of a shift in denudation rate from 𝜀1 to 𝜀2, but rather in terms of a change in the hillslope transport coefficient K, from K1 to K2. The values of K1 and K2 are defined for each cluster independently of the others. We perform a MCMC inversion to obtain the marginal posterior distributions for T and the various K, and assess the likelihood of the scenario (Figure 7). We have excluded basin 6 from the analysis, as it displays a clear outlier behavior in the relationship between denudation rate and hilltop curvature (Figure 3), which is likely due to a transition from transport- to weathering-limited conditions, as noted above.
Scenario 1: results of an inversion of the 10Be and 14Cin situ data assuming a step change in the hillslope transport coefficient K (from K1 to K2). Each cluster has different values of K1 and K2, but the timing of the change is synchronous. The different clusters, from 1 to 3, are identified by the colors blue, red and yellow, respectively. (A) Cumulative probability density curves for K1 and K2, for each cluster, before (solid line) and after (dashed line) the transition. The black line corresponds to the compilation of K values from Richardson et al. (2019). (B) Joint posterior distribution of K1 and K2. The black contours correspond to the probability density of the compilation of K values from Richardson et al. (2019). The thick solid grey line corresponds to K1 = K2. Thin grey solid and dashed lines correspond to various values of the ratio K2/K1, and the corresponding increase in denudation rate, from ×2 to ×100. (C) The orange distribution corresponds to the transition time between K1 and K2, which is common for all basins. The blue curve is the Greenland oxygen-isotope ice record of Grootes et al. (1993), shown for reference. (D) Comparison of observed and predicted 10Be and 14Cin situ concentrations. The black line indicates equality between the observed and predicted concentrations, and grey bands indicate deviations by ±20% and ±50%.
Scenario 1 implies an order-of-magnitude variation in K and denudation rates for clusters 1 and 2, and a more subdued evolution for cluster 3 (Figure 7A,B). Such large-amplitude changes in denudation rates have been reported by Hippe, Jansen, et al. (2021); however, in our modelling framework, they imply very high K2 values, outside the range reported by Richardson et al. (2019). This shift would have occurred after 30 ka BP, and most likely between 15 and 5 ka BP (orange colored distribution in Figure 7C). The inverted parameters predict most of the concentrations to within 20% of the observed values for both nuclides, with up to 6 catchments deviating more, but still within 50% of the observed values (Figure 7D). The very high values of K2 for two of the clusters imply that a variation in denudation through climatically-driven modulation of the diffusive hillslope transport coefficient K is not likely to fully explain the observed disequilibrium, and suggest that mass-wasting processes might be involved. However, we did not observe widespread recent landsliding in the studied basins.
In Scenario 2, we explore the instantaneous removal of a given thickness of soil H occurring at a given time T (Figure 8). Background denudation rates are set by using a constant value of the hillslope transport coefficients K, which are the same before and after the perturbation, and applying Equation (3). The values of H and K are defined for each cluster independently of the others. This inversion yields transport coefficients K which are within the range of the compilation by Richardson et al. (ibid.). The soil removal event is predicted to have occurred very recently, most likely during the last centuries (Figure 8C). The amplitude of removal is variable across clusters, with a depth up to 3 m at cluster 2, but <1 m for the others (Figure 8B). The removal of several meters of material may not be realistic, but it should be noted that the calculation uses a regolith density of 1.3 g/cm3, which may be an underestimate depending on the nature of the cover and the type of process involved. Scenario 2 is also successful at reproducing the observed concentrations (Figure 8D).
Scenario 2: results of an inversion of the 10Be and 14Cin situ data based on the scenario considering instantaneous soil removal at a given time. Each cluster has a different value of the thickness of soil removed. Each cluster has a different value of the hillslope transport coefficient K, which is constant through time and does not change at the transition in this scenario. The different clusters, from 1 to 3, are identified by the colors blue, red and yellow, respectively. (A) Cumulative probability density curves for K, for each cluster. The black line corresponds to the compilation of K values from Richardson et al. (2019). (B) Distribution of removed soil thickness for each cluster. (C) Distribution of the joint soil removal time. The blue curve is the Northern Hemisphere temperature anomaly from Briffa (2000), for reference. (D) Comparison of observed and predicted 10Be and 14Cin situ concentrations. The black line indicates equality between the observed and predicted concentrations, and grey bands indicate deviations by ±20% and ±50%.
5. Discussion
5.1. Multi-nuclide dataset analysis approach
Our reconstruction of the denudation history from the 14Cin situ and 10Be concentrations is inspired by the scenarios proposed by Hippe, Jansen, et al. (2021), with the difference that we incorporate information about the catchment morphology in our inversion, by using measured hilltop curvature and inverting for changes in the hillslope transport coefficient. We follow the method proposed by Godard, Salgado, et al. (2021), Godard, Siame, et al. (2024), Godard, Mudd, et al. (2025), which consists in formulating the time variations of surface processes in terms of efficiency parameters instead of individual denudation histories for each basin or site. This approach relies on the idea that denudation variations are strongly correlated, both in their timing, as was considered by Hippe, Jansen, et al. (2021), and in their amplitude, so that, for example, a twofold increase in the denudation rate of one catchment is likely to have affected nearby catchments with similar relative changes in intensity. Our approach uses a linear relationship between denudation rates and a morphological characteristic of the basin to decrease the dimension of the inversion problem. In our case we use Equation (3) with a hillslope transport coefficient considered common to all basins belonging to the same cluster. We could have used other geomorphological laws that express denudation rates as a function of a quantifiable morphological variable and an efficiency parameter setting the intensity of the considered process. In addition, our approach allows us to compare the inversion results with published values of K (Richardson et al., 2019).
Scenarios 1 and 2 consider one single event of either denudation acceleration or soil removal. Several sites in the FMC suggest various stages of surface processes as a response to climatic or anthropic perturbations, in particular over the Holocene (e.g. Surmely et al., 2009; Chassiot, Miras, et al., 2018; Chassiot, Simonneau, et al., 2022). Assuming a single event forces the inversion to maximize the amplitude of the surface process response, and the summation of a sequence of lower intensity events through time, either denudation acceleration or soil removal, might yield similar TCN signals. Such a sequence may be more likely than a single very intense event, but constraining its detailed timing would dramatically increase the number of parameters in our inversion, to the point of rendering it intractable. Additionally, the two scenarios we consider are not mutually exclusive, and we could speculate that the measured TCN signal results both from a change in denudation rates at the Pleistocene–Holocene climatic transition and from a recent, anthropogenically driven soil-removal event.
Our approach directly relies on the conditions of applicability of Equation (3), which include continuous hillslope transport. Consequently, Equation (3) cannot account for processes such as mass wasting and the discontinuous denudation of hillslope material, which may occur during drastic climatic changes. This might be a reason for the abnormally high recent transport coefficient (K2) obtained through the denudation change scenario (Figure 7), when compared with the range of values reported by Richardson et al. (2019).
We considered basins distributed across a geomorphic gradient with a strong contrast in relief and soil cover, but no differential in tectonic rock uplift. The k-mean clustering allowed us to identify subsets of catchments with comparable characteristics. This choice was guided by the diversity of observed 14Cin situ/10Be ratios in our dataset suggesting contrasted dynamics and responses to environmental changes. Basins with relatively low hillslope angles (<20°, clusters 1 and 3) clearly display the lowest amount of deviations between the two nuclides, with some catchments almost at equilibrium (Figures 1 and 2). The interior basins (cluster 1) show very limited impacts of climatically driven changes in hillslope dynamics, with almost no K variation in Scenario 1 (denudation-rate change; Figure 7) and very limited soil stripping in Scenario 2 (Figure 8). In contrast, cluster 2 in the Cévennes displays the largest variations, highlighting the sensitivity of this steep landscape to perturbations. This observed spatial gradient in sensitivity contrasts with some of the modeling results of Mudd (2016), who showed that, on the contrary, slowly eroding landscapes should display the highest response. We note that the catchments along the margin are also subject to a drastically different hydro-climatic regime compared with the interior of the massif, with occurrences of extreme rainfall events and large discharge variability (Desormeaux et al., 2022). The contrast in sensitivity between the two domains is then not only controlled by the intrinsic dynamics of the 14Cin situ–10Be system, but also by changes in the nature of surface processes through time with a possible intensification of extreme precipitation and discharge events along the topographic margin associated with the onset of a wetter climate during the Holocene (e.g. Peyron et al., 1998; Ponel et al., 2022).
5.2. Comparison with other 14Cin situ datasets
Only a limited number of previous studies using the 14Cin situ–10Be system in detrital sediments have been published, providing a glimpse into changes of Earth surface processes over the Late Quaternary (e.g. Kober et al., 2019; Hippe, 2017). Due to this limited number of previous studies and the geomorphic and climatic diversity of the settings they investigated, there is no single analysis framework for such datasets, which complicates the generalization of the results. We also note that these studies almost systematically observe significant deviations from steady-state ratios, which probably reflects the major variations in environmental parameters over the last tens of millennia (Godard, Mudd, et al., 2025), even though the actual processes invoked in each case are very different in nature and timing. In some cases, this combination of processes make the signal very difficult to interpret (Schmidt et al., 2026).
Hippe, Jansen, et al. (2021) showed a major disequilibrium of the 14Cin situ–10Be system in catchments of the Altiplano, which could be explained either by punctual events of denudation acceleration or soil removal. All these events occurred after 5 ka BP and Hippe, Jansen, et al. (ibid.) suggest that they could be related to wetter conditions over the area and coeval changes in agricultural practices, highlighting the connections between climatic variability and human activities. Our study converges towards similar interpretations but with different timings, as we observe a variable degree of disequilibrium in our basins that seems to be driven by how basins respond to perturbations depending on their hillslope gradient. For example, we could postulate that steep basins of the margin are close to threshold hillslope conditions and will display a shorter and non-linear response to changes in transport coefficient, whereas basins with lower slopes are more likely to stay in the linear domain, with more subdued responses (Godard and Tucker, 2021).
Slosson et al. (2022) observed similar low ratios and a strong disequilibrium of the 10Be–14Cin situ system in large basins draining the eastern flank of the Argentinian Andes, which they interpret as the manifestation of cumulated storage events on hillslopes within these basins. However, the setting of Slosson et al. (ibid.) is very different from the area we investigate: our basins are smaller, the glacial imprint is weak if not absent, and there are no active tectonics. We cannot directly compare their results with ours, as we have only a limited number of plausible processes, and transient storage is not realistic in our setting.
Fülöp, Codilean, et al. (2020) also obtained detrital 10Be–14Cin situ data in SE Australia, on basins much larger than the ones investigated in our study. They observe low ratios of variable amplitude which in their setting indicate complex sediment reworking in large-scale fluvial systems. Such processes are not relevant for the type of landscape and basins of our study. A recent study by Towers et al. (2025) focuses exclusively on 14Cin situ to investigate denudation processes in the postglacial landscape of the Scottish Highlands. They infer a significant contribution of paraglacial deposits in the studied catchment to the present-day sediment flux in the main channel, and thus highlight the lasting impact of glaciation on these processes. Because of a different setting without major glacial impact and different methods, we cannot conduct a meaningful comparison of their results with ours, but their study underscores the interest of 14Cin situ as a tracer of Holocene processes.
5.3. Comparison with existing constraints on recent climatic and anthropic forcings in the FMC
Although we selected our study area and the sampled basins to minimize the direct impact of glacial shielding and erosion during the Late Pleistocene, periglacial conditions were present in the landscape, with localized gelifraction and associated slope deposits. In many places, the development of periglacial formations on hillslopes was actually a key factor of sensitivity to later erosion events triggered by Holocene climatic fluctuations or anthropogenic land occupation changes (Dendievel et al., 2015). An extensive cover of clasts recently produced by gelifraction could be a factor promoting intense regolith removal during climatic changes, and could correspond to the major increase in denudation in Scenario 1.
Schaller, von Blanckenburg, et al. (2002) provided the only 10Be-based paleo-denudation dataset in the Massif Central, from terraces of the Allier and Dore rivers, located northwest of our study area. While their 30-ka record presents some minor fluctuations through time, they do not observe major changes in denudation rates between glacial and inter-glacial conditions. In particular, all their paleo-denudation rates are within the range of observed modern detrital 10Be denudation rates at similar locations (ibid.). Crucially, such sequences based solely on 10Be often display limited sensitivity to short period variations, in particular in slowly erosive settings (Schaller and Ehlers, 2006), and thus their results are not incompatible with ours, as we take advantage of short-lived nuclides such as 14Cin situ to detect fluctuations of potentially larger amplitude.
Human occupation has also been recognised, with scattered settlements and a variable imprint on vegetation throughout most of the Holocene (Cubizolle et al., 2014), although anthropogenic pressure increased markedly during the historical period (Surmely et al., 2009; Lavrieux et al., 2013; Jacob-Rousseau and Astrade, 2014). Notably, the eastern part of the FMC was largely deforested in the XVIIIth and XIXth centuries, facilitating soil erosion (Astrade et al., 2011). The extent of this deforestation and the recognition of the associated geomorphological crisis triggered a massive restoration program in French mountainous areas in the second half of the XIXth century and into the XXth century. A recent, historical, major increase in sediment fluxes associated with human-induced changes in vegetation cover is also well documented in lakes of the FMC (Macaire et al., 1997; Ballut et al., 2008). The recent erosion history of the eastern FMC seems to be influenced by the general patterns of climatic evolution over western Europe, with episodic interference from human activities. As proposed by Defive (2013), both forcings act in conjunction, with their timing mostly dictated by climate but their spatial expression controlled by human land use.
The two types of forcings could correspond to the end-member scenarios we consider in our inversion of the dataset. In the hypothesis of a single increase in denudation rates (Scenario 1), the timing of the associated event appears to be restricted to the 15–5 ka BP interval. Within this time frame, one candidate is the Pleistocene–Holocene transition, whose impact on the acceleration of TCN-derived denudation rates worldwide has recently been proposed by Godard, Mudd, et al. (2025). Some lacustrine records of the eastern part of the FMC extend into the Pleistocene and display evidence for an increase in sedimentation rates around 10 ka (Martin et al., 2019). However, the limited number of sites prevents us from concluding whether this is a regional pattern, and changes in the nature of the sediments make the interpretation of this signal ambiguous. Other records suggest instead a stabilizing influence of the Holocene expansion of vegetation (Degeai and Pastre, 2009) or a decrease in sedimentation rates (Thouveny et al., 1994). At other sites in the Ardèche mountains, pulses of non-organic matter in the sediment record are unambiguously associated with climatic perturbation phases during the Late Pleistocene and Early Holocene (Dendievel et al., 2015). While an increase in denudation rates (Scenario 1) around this transition could explain our observed concentrations, the associated increase in the transport coefficient is too high to be compatible with our assumption of continuous hillslope diffusive transport. The response we observe cannot derive solely from changes in this transport coefficient, as has been proposed for some systems (Hughes et al., 2009), and we expect a contribution from discontinuous, stochastic mass wasting.
The FMC sites that provide various records of Holocene environmental changes suggest variable erosion impacts from anthropic pressure, with changes connected to climatic and socio-economic factors (e.g. Chassiot, Simonneau, et al., 2022; Lavrieux et al., 2013). The two forcings are highly interconnected, and their respective contributions are usually difficult to deconvolve, hindering identification of the impact of individual climatic events such as the 8.2 ka cold and wet event (Magny et al., 2003). Climatic and anthropogenic perturbations can contribute to the exhumation of deep soil and bedrock and deliver material with lower 14Cin situ/10Be ratios in the fluvial system. We note, however, that the soil thicknesses required to explain these ratios are >2 m in one of the clusters, corresponding to the Cévennes area. Although deforestation was intense in parts of this landscape, such a substantial soil removal depth is clearly unrealistic. Conversely, this scenario is permissible for the other clusters, with soil-removal thicknesses <1 m. In the end, both end-member scenarios are probably too extreme to independently explain the depressed ratios we observe, and we suggest a combination of different processes, such as the activation of hillslope processes around the Pleistocene–Holocene transition and the recent anthropogenic disturbances related to the XVIIIth–XIXth deforestation.
6. Conclusion
We presented a new dataset of 14Cin situ concentrations in river sediments complementing existing 10Be data for 10 basins along the eastern margin of the French Massif Central. We observe varying degrees of disequilibrium between the two nuclides, which point to a transient evolution of surface processes over the Holocene or the Late Glacial period. The diversity in 14Cin situ/10Be ratios suggests complex spatial variations in the nature and timing of these changes, and different responses across our study area, in part related to the contrasting geomorphological contexts of the steep margin and the low-relief interior of the Massif. We explored various scenarios to reproduce the observed concentrations, including a recent increase in denudation rates and a soil-removal event.
Among previous studies that have highlighted a disequilibrium in the 10Be–14Cin situ system for river sediments, our work is the first to incorporate topographic information into the inversion of TCN concentrations. This approach accounts for similarities in denudation dynamics between adjacent basins while reducing the number of inferred parameters. A scenario implying an increase in denudation rates around 10 ka can explain most of the observed concentrations, but would suggest hillslope transport coefficients that are outside the range of the values reported in the literature. Alternatively, a very recent (within the last few centuries) removal of soil or regolith could also explain the observed concentrations and disequilibrium, with variable depths of removal across the studied area. A combination of both processes—variations in background denudation and soil removal—most likely involving multiple events, is probably driving the transient evolution of surface processes in this area.
14Cin situ is now an analytically mature TCN system with robust systematics, even if the low number of extraction facilities limits its broader use. Interpreting the widely observed disequilibrium of the 10Be–14Cin situ remains very challenging and context-specific, as the system is highly sensitive to recent changes in surface processes, even over very short timescales of a few centuries. In our study and others, such as Hippe, Jansen, et al. (2021), a panel of different scenarios can be invoked to explain the observed deviations in isotopic ratios. Addressing these ambiguities requires methodological developments in interpretating this signal and also integration of other constraints on the recent evolution of surface processes in a unified analysis framework.
Acknowledgments
This research was supported by ANR TOPOEXTREME (ANR-18-CE01-0017), INSU/Tellus (Maraca) and Radiate (20002215-ST-1.1) grants. We thank two anonymous reviewers for constructive and insightful comments.
Declaration of interests
The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.