# Determining the Plio-Quaternary uplift of the southern French Massif Central; a new insight for intraplate orogen dynamics

Oswald Malcles, Philippe Vernant, Jean Chery, Pierre Camps, Gaël Cazes, Jean-François Ritz, David Fink

► To cite this version:

Oswald Malcles, Philippe Vernant, Jean Chery, Pierre Camps, Gaël Cazes, et al.. Determining the Plio-Quaternary uplift of the southern French Massif Central; a new insight for intraplate orogen dynamics. Solid Earth, 2020, 11 (1), pp.241-258. 10.5194/se-11-241-2020. hal-04749762

HAL Id: hal-04749762

https://hal.science/hal-04749762v1

Submitted on 23 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Determining the Plio-Quaternary uplift of the southern French Massif Central; a new insight for intraplate orogen dynamics

Oswald Malcles<fnref n="1" />, Philippe Vernant<fnref n="1" />, Jean Chéry<fnref n="1" />, Pierre Camps<fnref n="1" />, Gaël Cazes<fnref n="2" />,<fnref n="3" />, Jean-François Ritz<fnref n="1" />, and David Fink<fnref n="3" />

**Correspondence:** Oswald Malcles (oswald.malcles@umontpellier.fr)

Received: 29 May 2019 – Discussion started: 11 June 2019

Revised: 16 January 2020 – Accepted: 21 January 2020 – Published: 26 February 2020

**Abstract.** The evolution of intraplate orogens is still poorly understood. Yet, it is of major importance for understanding the Earth and plate dynamics, as well as the link between surface and deep geodynamic processes. The French Massif Central is an intraplate orogen with a mean elevation of 1000 m, with the highest peak elevations ranging from 1500 to 1885 m. However, active deformation of the region is still debated due to scarce evidence either from geomorphological or geodetic and seismologic data. We focus our study on the southern part of the Massif Central, known as the Cévennes and Grands Causses, which is a key area to study the relationship between the recent geological deformation and landscape evolution. This can be done through the study of numerous karst systems with trapped sediments combined with the analysis of a high-resolution digital elevation model (DEM).

Using the ability of karst to durably record morphological evolution, we first quantify the incision rates. We then investigate tilting of geomorphological benchmarks by means of a high-resolution DEM. We finally use the newly quantified incision rates to constrain numerical models and compare the results with the geomorphometric study.

We show that absolute burial age ($^{10}$Be/$^{26}$Al on quartz cobbles) and the paleomagnetic analysis of karstic clay deposits for multiple cave system over a large elevation range correlate consistently. This correlation indicates a regional incision rate of $83 \pm 17/\sqrt{5}$ m Ma$^{-1}$ during the last ca. 4 Myr (Pliocene–Quaternary). Moreover, we point out through the analysis of 55 morphological benchmarks that the studied region has undergone a regional southward tilting. This tilting

is expected as being due to a differential vertical motion between the northern and southern part of the studied area.

Numerical models show that erosion-induced isostatic rebound can explain up to two-thirds of the regional uplift deduced from the geochronological results and are consistent with the southward tilting derived from morphological analysis. We presume that the remaining unexplained uplift is related to dynamic topography or thermal isostasy due to the Massif Central Pliocene–Quaternary magmatism.

Integrating both geochronology and morphometrical results into lithospheric-scale numerical models allows a better understanding of this intraplate–orogen evolution and dynamic. We assume that the main conclusions are true to the general case of intraplate deformation. That is to say, once the topography has been generated by a triggering process, rock uplift is then enhanced by erosion and isostatic adjustment leading to a significant accumulation of mainly vertical deformation.

## 1 Introduction

For the past few decades, plate-boundary dynamics have been, to a first order, well understood. This is not the case for intraplate regions, where short-term ($10^3–10^5$ years) regional strain rates are low and the responsible dynamical processes are still in debate (e.g., Calais et al., 2010, 2016; Vernant et al., 2013; Tarayoun et al., 2017). Intraplate deformations evidenced by seismic activity are sometimes explained by a transient phenomenon (e.g., glacial isostatic rebound or

<footnote n="1">¹Geosciences Montpellier, CNRS and University of Montpellier, Montpellier, France</footnote>

<footnote n="2">²School of Earth and Environmental Sciences, University of Wollongong, Wollongong, Australia</footnote>

<footnote n="3">³Australian Nuclear Science and Technology Organisation, Lucas Heights, Australia</footnote>
<figure><img src="image_3.png" /></figure>

**Figure 1.** A 30 m resolution DEM of the French Massif Central and grey scale underlying the DEM is the DEM derived slope map. Examples of finite incision typical of the French Massif Central in (a) crystalline basement (Seuge canyon) and (b) limestone plateau (Tarn canyon). Location of the study area in red box (Fig. 9) and enumerated site (1) is the Rieutord canyon (43.958° N, 3.709° E) where TCN measurements were taken; (2) and (3) are the Leicasse cave system (43.819° N, 3.56° E), and the Garrel cave system (43.835° N, 3.616° E), respectively, where paleomagnetic analyses have been done, and (4) is the Lodève basin (43.669° N, 3.382° E) with dated basaltic flows. Bottom panel is an example of typical topographic profile used for the numerical model set up. Note the southwestern area with large plateau dissected by a canyon, and the rugged area with a steep valley is called the Cévennes. They are typical regional limestone and crystalline morphology, respectively.
hydrological loading). However, to explain the persistence through time of intraplate deformation and explain the high finite deformation that we can observe in the topography in many parts of the world as for instance the Ural Mountains in Russia, the Blue Mountains in Australia, or the French Massif Central, one needs to invoke continuous processes at the geological timescale.

Located in the southwestern Eurasian plate (Fig. 1), the French Massif Central is an ideal case to study these processes because a high-resolution digital elevation model (DEM) encompasses the whole region, and widespread karstic areas are present along its southern and western edges, allowing the possibility to quantification of the landscape evolution rates thanks to terrestrial cosmogenic nuclides (TCN) burial ages. The region is characterized by a mean elevation of 1000 m with summits higher than 1500 m. Such a topography is likely to be the result of recent active uplift, and as the Cévennes mountains experiences an exceptionally high mean annual rainfall (the highest peak, Mont Aigoual, records the highest mean annual rainfall in France of 4015 mm) it raises the question of a possible link between erosion and uplift as previously proposed for the Alps (Champagnac et al., 2007; Vernant et al., 2013; Nocquet et al., 2016). This region currently undergoes a small but discernible deformation, but no significant quantification can be deduced due to the scarcity in seismicity (Manchuel et al., 2018). In addition, GPS velocities are below the uncertainty threshold of GPS analyses (Nocquet and Calais, 2003; Nguyen et al., 2016).

In this study we focus on the Cévennes mountains and the Grands Causses (Fig. 1) area, where cave systems with trapped sediments are known over a widespread altitude range. South and west of the crystalline Cévennes mountains, prominent limestone plateaus, named Grands Causses, rise to 1000 m and are dissected by few canyons that are several hundreds of meter deep. The initiation of incision, its duration, and the geomorphic processes leading to the present-day landscape remain poorly constrained. A better understanding of the processes responsible for this singular landscape would bring valuable information on intraplate dynamics, especially where large relief exists.

## 2 Geological background

The oldest rock units in the study area were formed during the Variscan orogeny (late Paleozoic, ~ 300 Ma; Brichau et al., 2007) and constitute the crystalline basement of the Cévennes. Between 200 and 40 Ma (Mesozoic and middle Cenozoic), the region was mainly covered by the sea, ensuring the development of an important detrital and carbonate sedimentary cover, which can reach several kilometers of thickness in some locations (Sanchis and Séranne, 2000; Barband et al., 2001). During the Mesozoic Era, an episode of regional uplift and subsequent erosion and alteration (called

the Durancian event) is proposed for the origin of the flat, highly elevated surface that persists today across the landscape (Bruxelles, 2001; Husson, 2014).

The area is affected by the major NE–SW trending Cévennes Fault system, a lithospheric-scale fault, inherited from the Variscan orogen. This fault system was reactivated several times (e.g., as a strike-slip fault during the Pyrenean orogen or as a normal fault during the Oligocene extension). During the Pyrenean orogeny, between 85 to 25 Ma (Tricart, 1984; Sibuet et al., 2004), several fault and fold systems affected the geological formations south of the Cévennes Fault, while very few deformations occurred farther north within the Cévennes and Grands Causses areas (Arthaud and Laurent, 1995). Finally, the Oligocene extension (~ 30 Ma) led to the counterclockwise rotation of the Corso-Sardinian block and the opening of the Gulf of Lion, reactivating some of the older compressive structures as normal faults. The main drainage divide between the Atlantic Ocean and the Mediterranean Sea is located in our study area and is inherited from this extensional episode (Séranne et al., 1995; Sanchis and Séranne, 2000).

Afterwards, during the Pliocene–Quaternary period, intense volcanic activity affected the region, from the Massif Central to the Mediterranean shoreline. This activity is characterized by several volcanic events that are well constrained in age (Dautria et al., 2010). The last eruption occurred in the Chaîne des Puys during the Holocene (i.e., the past 10 kyr; Nehlig et al., 2003; Miallier et al., 2004). Some authors proposed that this activity was related to a hotspot underneath the Massif Central that led to an observed positive heat-flow anomaly and a possible regional Pliocene–Quaternary uplift (Granet et al., 1995; Barruol and Granet, 2002). Geological mapping at a different scale can be found at: http://infoterre.brgm.fr/ (last access: 21 February 2020).

Despite this well-described overall geological evolution, the onset of active incision that has shaped the deep valleys and canyons (e.g., Tarn or Vis river, Fig. 1) across the plateaus and the mechanisms that controlled this incision are still debated. One hypothesis proposes that canyon formation was driven by the Messinian salinity crisis with a drop of more than 1000 m in Mediterranean Sea level (Moccochain, 2007). This, however, would not explain the fact that the Atlantic watersheds show similar incision. Other studies have suggested that the incision is controlled by the collapse of cave galleries that lead to fast canyon formation mostly during the late Quaternary, thus placing the onset of canyon formation at only a few hundreds of thousands of years ago (Corbel, 1954). More recently, it has been proposed (based on relative dating techniques and sedimentary evidence) that incision during the Quaternary was negligible (i.e., less than a few tens of meters) and that the regional morphological structures seen today occurred around 10 Ma (Séranne et al., 2002; Camus, 2003).
## 3 Materials and methods

In this paper, we provide new quantitative constraints on both the timing of incision and the rate of river downcutting in the central part of the Cévennes and the Grands Causses that has resulted in the large relief between plateau and channel bed.

We employ two methods to infer allochthonous karstic infilling age and associated river downcutting. First, we use quartz cobbles to measure the concentration of cosmogenic $^{10}$Be and $^{26}$Al isotopes. The $^{10}$Be/$^{26}$Al ratios provide burial ages of these karstic infillings. Second, paleomagnetic analyses of clay deposits are used to obtain paleo-polarities. In both cases, vertical profiles among tiered caves systems and horizontal galleries could provide local incision rate information. By analyzing a high-resolution DEM (5 m), we show that the region is affected by a southeastward regional tilting. Our results allow quantification of the role of the Pliocene–Quaternary incision on the Cévennes landscape evolution and constraint of the numerical modeling from which we derive the regional uplift rates and tilt of geomorphological markers.

If incision is initiated by uplift centered on the north of the area where elevations are at a maximum, it will lead to tilting of fossilized topographic markers as strath terraces. Our research approach provides an opportunity to discriminate between three possible explanations for the current terrain morphology. The first is based on old uplift and old incision (Fig. 2a). In this case, apparent incision rates would be very low. For instance, if incision commenced 10 Ma (Sérane et al., 2002), we would find surface tilting, but cosmogenic burial dating with $^{10}$Be/$^{26}$Al, which cannot discern ages older than ~5 Ma due to excessive decay of $^{26}$Al, would not be possible. The second possibility (Fig. 2b) is that the uplift is old, and incision consequently follows but with a time lag. Here the incision rate would be rather fast, but no tilting is expected for the river-related markers because no differential uplift occurs after their formation. Finally, the third possibility (Fig. 2c) is that uplift and incision are concurrent and recent (i.e., within the timescale of cosmogenic burial dating), and thus we would expect burial ages < 5 Myr, relatively high incision rates, and tilting of morphological markers. These different proposals for the temporal evolution of the region will be compared using numerical modeling.

# 4 Determining the incision rates in the Cévennes and the Grands Causses regions

## 4.1 Principles and methods

### 4.1.1 Karst model

No evidence of important aggradation events has been reported in the literature for the studied area. Therefore, we base our analysis on a per descensum infill model of the

karst networks whereby sediments are transported and then deposited within cave galleries close to base level. When cave systems and entry passages are near the contemporaneous river channel elevation (including higher levels during floods), the deposition into caves of sediments, from clay to cobbles, occurs, especially during flood events. Subsequent river incision into bedrock creates a relative base level drop (due to uplift or sea level variations). The galleries associated with the former base level are now elevated above the new river course and become disconnected from further deposition. Hence fossilized and trapped sediments throughout the cave network represent the cumulative result of incision. In this commonly used model (Granger et al., 1997; Audra et al., 2001; Stock et al., 2005; Harmand et al., 2017), the higher the gallery elevation (relative to the present-day base level), the older the deposits in that gallery. As a result, the objective here is to quantify a relative lowering of the base level in the karst systems, with the sediments closest to the base level being the youngest deposits, and note that we do not date the cave network creation which may very well precede river sediment deposition.

Within individual canyons, successions of gallery networks across the full elevation range, from plateau top to modern river channel, were not always present, and often sampling could not be conducted in a single vertical transect. Thus, we make the assumption of lateral altitudinal continuity, i.e., that within a watershed, which may contain a number of canyons, the sediments found in galleries at the same elevation were deposited at the same time. Inside one gallery, we use the classical principle of stratigraphy sequence (i.e., the older deposits are below the younger ones). More information and detailed relationships concerning the karstic development and geometric relationship between karstic network and morphological markers can be found in Camus (2003). In any cases, our aim is not to date the formation of the galleries or to explain the formation processes (e.g., past preferential alteration layer) but to use the time information brought by the sediments that have been trapped in the cave system. Therefore, we apply the commonly used model (example in Harmand et al., 2017) that has been proven by Granger et al. (1997, 2001). For the cave topographic survey, we refer the reader to https://data.oreme.org/karst3d/karst3d_map (last access: 21 February 2020), which provides a 3-D survey.

### 4.1.2 Burial ages

Burial dating using TCN is now a common tool to quantify incision rates in a karstic environment (Granger and Muzikar, 2001; Stock et al., 2005; Moccochain, 2007; Tassy et al., 2013; Granger et al., 2015; Calvet et al., 2015; Genti, 2015; Olivetti et al., 2016; Harmand et al., 2017; Rovey II et al., 2017; Rolland et al., 2017; Sartégou, 2017; Sartégou et al., 2018a). This method relies on the differential decay of TCN in detrital rocks that were previously exposed to cosmic ra-
<figure><img src="image_4.png" /></figure>

**Figure 2.** Conceptual models for landscape evolution. Top panel is the initial stage (prior to uplift). Each panel represents a possible scenario explaining current morphology: (a) old uplift and old incision, (b) old uplift and recent incision, and (c) both recent uplift and incision. Blue arrows and associated ages show expected result (or absence of) for burial dating. Red level represents morphological markers that are fossilized when reaching the surface, accumulating afterward (or not) the differential uplift by finite tilting.

diation before being trapped in the cave system. With this in mind, the ¹⁰Be and ²⁶Al nuclide pair is classically used as (i) both nuclides are produced in the same mineral (i.e., quartz), (ii) their relative production ratio is relatively well constrained (we use a standard ²⁶Al/¹⁰Be pre-burial ratio of 6.75; see Balco et al., 2008), and (iii) their respective half-lives (about 1.39 and 0.70 Myr for ¹⁰Be and ²⁶Al, respectively) are well suited to karstic and landscape evolution study, with a useful time range of ~100 kyr to ~5 Myr.

To quantify the incision rate of the limestone plateau of the Cévennes area, we analyzed quartz cobble infillings from four caves of the Rieutord canyon (Fig. 1), this canyon is well suited for such a study because horizontal cave levels are tiers over 200 m above the current river level and are directly connected to the canyon, leading to a straight relationship between river elevation and the four cave infillings that we have sampled (Cuillère Cave, Route Cave, Camp de Guerre Cave and Dugou Cave). Furthermore, the source of the cobbles is well known and identified as the upstream part of the Rieutord river, some tens of kilometers northward, providing a unique sediment origin composed of granite and metamorphic rocks embedded with numerous quartz veins. All samples (Example Fig. 3) were collected far enough away (> 20 m) from the cave entrance and deep enough below the surface (> 30 m) to avoid secondary in situ cosmogenic production of ¹⁰Be and ²⁶Al in the buried sediments.

The quartz cobbles were first crushed and purified for their quartz fraction by means of sequential acid attack with aqua regia (HNO₃ + 3HCl) and diluted hydrofluoric acid (HF).

<figure><img src="image_5.png" /></figure>

**Figure 3.** Example of quartz cobbles sampled for burial dating. Location: Cuillère Cave (Site 1, Fig. 1).

Samples were then prepared according to ANSTO's protocol (see Child et al., 2000), and ~300 µg of a ⁹Be carrier solution was added to the purified quartz powder before total dissolution. Accelerator mass spectrometry (AMS) measurements were performed on the 6MV SIRIUS AMS instrument at ANSTO, and results were normalized to KN−5−4 (for Be,
see Nishiizumi et al., 2007) and KN-4-2 (for Al) standards. Uncertainties for the final $^{10}$Be and $^{26}$Al concentrations include AMS statistics, 2% (Be) and 3% (Al) standard reproducibility, 1% uncertainty in the Be carrier solution concentration, and 4% uncertainty in the natural Al measurement made by inductively coupled plasma optical emission spectrometry (ICP-OES), in quadrature. Sample-specific details and results are found in Table 1.

### 4.1.3 Paleomagnetic analysis

In parallel with burial dating, we analyzed the paleomagnetic polarities within endokarstic clay deposits within two main cave systems: the Grotte-exsurgence du Garrel (Garrel) and the Aven de la Leicasse (Leicasse) (Fig. 1). These two cave systems allowed us to collect samples along a more continuous range of elevations than the one provided by the Rieutord samples (for burial age determination) and also to extend the spatial coverage to the southern Grands Causses region. Thanks to the geometry of these two cave systems, we sampled a 400 m downward base level variation. The sampling was done along vertical profiles from a few tens of centimeters to 2 m high by means of Plexiglas cubes with a 2 cm edge length (Fig. 4), which are used as a pastry cutter. We were not able to analyse clay samples from Rieutord canyon because no reliable clay infilling was found in the Rieutord caves.

Demagnetization was performed with an applied alternative field up to 150 mT using a 2G760 cryogenic magnetometer, equipped with the 2G600 degaussing system controller. Before this analysis, each sample remained at least 48 h in a null magnetic field, preventing a possible low-coercivity viscosity overprinting the detrital remanent magnetization (DRM) (Hill, 1999; Stock et al., 2005; Hajna et al., 2010). If the hypothesis of instantaneously locked-in DRM seems reasonable compared with the studied time span, it is important to keep in mind that the details of DRM processes (as for instance the locked in time) are not well understood (Tauxe et al., 2006; Spassov and Valet, 2012) and could possibly lead to small variations (few percents) in the following computed incision rates.

Because fine clay particles are expected to be easily reworked in the cave, careful attention was paid to the site selection, and current active galleries were avoided. Clays deposits had to show well-laminated and horizontal layering in order to prevent analysis of in situ-produced clay (from calcification) or downward drainage by an underneath diversion gallery that could strongly affect the obtained inclination (and also the declination to a minor extent). Note that for the paleo-polarity study alone, small inclination or declination variations will not result in false polarities.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Example of clay sampling for the paleomagnetic study. Location at the entrance shaft (highest elevation of every samples, ~ 580 m a.s.l.; Leicasse cave system, Site 2, Fig. 1).</figcaption></figure>

## 4.2 Quantifying the average incision rates

### 4.2.1 Local incision rate from burial ages (Rieutord canyon)

The relationship between burial ages and incision is shown in Fig. 5. For the four caves, we observed a good relationship between burial ages and finite incision, except for the Camp de Guerre Cave (CDG) site; the higher the cave is, the older the burial ages are. Burial ages for the Cuillère Cave, Dugou Cave, Camp de Guerre Cave and Route Cave are 2.16±0.15, 0.95±0.14, 0.63±0.1, and 0.21±0.1 Myr, respectively. This is consistent with the supposed cave evolution and first-order constant incision of the Rieutord canyon. The CDG age has to be considered with caution. The CDG cave entrance located in a usually dry thalweg can act as a sinkhole or an overflowing spring depending on the intensity of the rainfall. The sample was collected in a gallery that showed evidence of active flooding ~ 10 m above the Rieutord riverbed; therefore the older than expected age, given the elevation of the cave, is probably due to cobbles that came from upper galleries during flood events. Forcing the linear regression to go through the origin leads to an incision rate of 83±35 m Ma⁻¹.
<table><caption>**Table 1.** Samples analytical results and parameters. Cave abbreviations are RTE for the Route Cave, CDG for the Camp de Guerre Cave, DUG for the Dugou Cave and CUI for the Cuillère Cave. Main parameters are the geographical coordinate (lat and long in decimal degrees), the elevation (a.s.l.), and the height (a.b.l., computed relatively to the surface river elevation). The concentration (atoms g⁻¹ quartz) of ¹⁰Be and ²⁶Al in collected sand samples are all AMS ¹⁰Be/Be and ²⁶Al/Al isotopic ratios, corrected for full procedural chemistry blanks and normalized to KN-5-4 and KN-4-2, respectively. The error in the parenthesis is for total analytical error in final average ¹⁰Be and ²⁶Al concentrations based on statistical counting errors in final ¹⁰Be/Be (²⁶Al/Al) ratios measured by AMS in quadrature, with a 1% error in ⁹Be spike concentration (or a 4% error in ²⁷Al assay in quartz) and a 2% (or 3%) reproducibility error based on repeat of AMS standards. Burial age (minimum) assuming no post-burial production by muons at given depth (all deeper than 30 m) in cave below surface and assuming initial ²⁶Al/¹⁰Be ratio is given by the production ratio of 6.75. The burial age error determined by using a ±1σ range in the measured ²⁶Al/¹⁰Be ratio.</caption><thead><tr><th>Cave</th><th>Lat</th><th>Long</th><th>Elevation</th><th>Height</th><th>Conc ¹⁰Be<br>(10<sup>4</sup> atoms g<sup>-1</sup>)</th><th>σ ¹⁰Be<br>(10<sup>3</sup> atoms g<sup>-1</sup>)</th><th>Conc ²⁶Al<br>(10<sup>5</sup> atoms g<sup>-1</sup>)</th><th>σ ²⁶Al<br>(10<sup>4</sup> atoms g<sup>-1</sup>)</th><th>²⁶Al/¹⁰Be<br>(and σ)</th><th>Burial age<br>and σ<br>(mega-annums)</th></tr></thead><tbody><tr><td>RTE</td><td>43.960</td><td>3.707</td><td>175</td><td>8</td><td>3.54</td><td>1.18</td><td>2.16</td><td>1.47</td><td>6.11 ± 0.46</td><td>0.20 ± 0.15</td></tr><tr><td>CDG</td><td>43.955</td><td>3.710</td><td>185</td><td>10</td><td>8.87</td><td>3.12</td><td>4.29</td><td>3.28</td><td>4.83 ± 0.41</td><td>0.67 ± 0.16</td></tr><tr><td>DUG</td><td>43.957</td><td>3.711</td><td>245</td><td>115</td><td>1.27</td><td>5.68</td><td>0.529</td><td>0.636</td><td>4.15 ± 0.53</td><td>0.99 ± 0.25</td></tr><tr><td>CUI</td><td>43.959</td><td>3.711</td><td>354</td><td>175</td><td>1.70</td><td>7.14</td><td>0.375</td><td>0.528</td><td>2.20 ± 0.32</td><td>2.28 ± 0.28</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Correlation diagram of finite incision and burial age for the Rieutord canyon (Site 1, Fig. 1). Finite incision is the elevation of the sampling site relative to the current riverbed. RTE for Route Cave, CDG for Camp de Guerre Cave, DUG for Dugou Cave and CUI for Cuillère Cave.</figcaption></figure>

These results show that at least half of the 300 m deep Rieutord canyon is a Quaternary incision. Extrapolating the obtained rate yields an age of 4.4 ± 1.9 Ma for the beginning of the canyon incision, which suggests that the current landscape has been shaped during the Pliocene–Quaternary period. To extend our spatial coverage and bring stronger confidence in our results, we combine Rieutord burial ages with paleomagnetic data from watersheds located on the other side of the Hérault watershed.

### 4.2.2 Local incision rate from paleomagnetic data (southern Grands Causses)

A total of 100 clay-infilling samples distributed over 13 sites (i.e., profiles) was studied. The lowest sample elevation above sea level (a.s.l.) is in the Garrel cave system (ca. 190 m) and the highest is in the Leicasse cave system (ca. 580 m a.s.l.). In the Leicasse cave system, we sampled

eight profiles totaling 60 samples. Profiles elevations are located between ca. 200 and ca. 400 m above the base level (a.b.l.), which corresponds to the elevation of the Buèges river spring at 170 m a.s.l.

In the Garrel cave system, we sampled five profiles for a total of 40 samples that range between 20 and 80 m a.b.l., defined by the Garrel spring at 180 m a.s.l. Given the very marginal difference in elevation between the local base levels from these two caves, we assume that they have the same local base level. At each studied site, if all the profile samples have the same polarity, the site is granted with the same polarity, either normal or reverse. If not (i.e., the profile displays normal and reverse polarities), we consider it as a transitional site. Figure 6 shows the results plotted with respect to the paleomagnetic scale (x axis) for the past 7 Ma and their elevation above the base level (y axis). The measured paleomagnetic polarities on each site are plotted several times for given incision rates supposed to be constant through time (this al-
<figure><img src="image_8.png" /></figure>

**Figure 6.** Constraining the incision rate in the Cévennes margin, using paleomagnetic polarities from clay deposits (black, grey, and white symbols) and burial ages (red crosses). Circles are from the Leicasse cave system, with LGP being the les Gours sur Pattes profile (see text); squares are from the Garrel cave system. Black, grey, and white symbols correspond to normal, transitional, and reverse polarities, respectively. Straight black lines define possible incision rates that are supposed stable through time (numbers in white rectangles define the correlation factor (Cf) values between the measured paleomagnetic polarities and the predicted paleomagnetic scale (see also Fig. 8)). Green hexagons show the U–Th ages obtained on speleothems in the Garrel cave system by Camus (2003).

lows the determination of different age models and analysis of their correlation with the distribution of paleomagnetic data; see below). First, we note a good agreement between samples located at the same elevation and being part of the same stratigraphic layer (Camus, 2003). This syngenetic deposition allows, as best explanation to prevent from a possible partial endokarstic reworking. Second, the different elevations of the galleries where we collected the samples allow us to propose that the Leicasse deposits encompass at least three chrons, while the Garrel deposits encompass only one. Third, one sample presents a transitional signal, and this sample that is located (geometrically) between one with a reversal signal (lower samples) and one with a normal signal (upper ones) is observed at les Gours sur Pattes (LGP) sampling site (Fig. 7). This provides a strong constraint on the age of the sediment emplacement in the Leicasse cave system with respect to the magnetostratigraphic timescale (Fig. 6).

Compared to the Leicasse cave system, the elevation/polarity results for the Garrel cave system are less constrained. Only one site shows a reverse polarity at 90 m a.b.l., and the transitional polarity found at 40 m a.b.l. is unclear (see tables in the Supplement). The rest of the polarities (72 samples) are all normal. Given that a U–Th age younger than 90 kyr was obtained for two speleothems (Camus, 2003) covering our samples collected at 40 m a.b.l. (Fig. 6), we consider that the emplacement of the clay deposits occurred during the most recent normal period and are therefore younger than 0.78 Ma (Fig. 6). The transition between the highest nor-

mal sample and the reversed one is located somewhere between 78 and 93 m a.b.l., suggesting a maximum base level lowering rate of $109 \pm 9 \text{ m Ma}^{-1}$.

To go further in the interpretation of our data and better constraint the incision rate, we performed a correlation analysis between observed and modeled polarities for a 0 to $200 \text{ m Ma}^{-1}$ incision rate range (linear rate, each $1 \text{ m Ma}^{-1}$). Modeled polarities are found using the intersection between sample elevation and incision rate line.

We obtained 10 possible incision rates with the same best correlation factor (Fig. 8) spanning from 43 to $111 \text{ m Ma}^{-1}$ (mean of $87 \pm 24 \text{ m Ma}^{-1}$). Taking into account the transitional signal of the LGP site in the Leicasse cave system yields a linear incision rate of $83 \pm 17/ - 5 \text{ m Ma}^{-1}$. Proposed uncertainties are based on previous and next transition-related estimated incision rate.

Using a similar approach for the Rieutord crystalline samples, we computed for the same incision-rate space, the distance, in a least-squares sense, between the modeled ages and the measured ones in order to check the cost function shape and acuteness. With this method, we determined a linear incision rate of $85 \pm 11 \text{ m Ma}^{-1}$ (Fig. 8). Those two results, based on independent computations, suggest the same first-order incision rate for the last 4 Ma of $84 \pm 21/ - 12 \text{ m Ma}^{-1}$. Given that the Rieutord, Garrel and Buèges rivers are all tributaries of the Hérault River, we propose that this rate represents the incision rate for the Hérault River watershed, in-
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Zijderveld diagram for three samples from les Gours sur Pattes (Leicasse, Site 2, Fig. 1) site. Stratigraphical order is from (a) (the older, base of the profile) to (c) (the younger, top of the profile).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Best incision rates based on paleomagnetic data (blue) and burial ages (red). The blue curve is the normalized smoothed (10 m Ma⁻¹) sliding window for better visualization) correlation between theoretical and observed polarities. The highest correlation corresponds to the best incision rates. The red curve is the RMSE for the linear regression through the burial ages data set shown on Fig. 4.</figcaption></figure>
ducing approximately 300–350 m of finite incision over the Pliocene–Quaternary period.

If the landscape is at first order in an equilibrium state, that is to say, if we preclude our incision rates from being a regressive erosional signal, the incision needs to be balanced by an equivalent amount of uplift. If the uplift rate is roughly correlated to the regional topography, lowest uplift rates would be expected in the south of our sampling sites, inducing regional tilting of morphological benchmarks. In the next part, we search for such evidence, which would suggest a differential uplift.

# 5 Geomorphometric signature

## 5.1 Tested hypothesis and methods

According to the Massif Central centered uplift hypothesis, morphological markers such as strath terraces, fluvio-karstic surfaces, or abandoned meanders should display a southward tilting due to differential uplift between the northern and the southern part of the region.

To investigate these differential vertical movement signals, we used the morphological markers available in the study area (Fig. 9). We used a 5 m resolution DEM analysis to identify the markers corresponding to surfaces with slope < 2°. This cutoff slope angle prevents the identification of surfaces related to local deformation such as a landslide or sinkhole. We point out that surface slope increase through time (e.g., apparent tilting) could be due to diffusion processes and unrelated to differential vertical displacements. However, that problem is addressed by the automatic selection and correction and the final manual check for residue random distribution (see below). The local river slope is on the order of 0.1°, so the 2° cutoff angle is far from precluding the identification of tilted markers. We also use a criterion based on an altitudinal range for a surface. This altitudinal span is set individually for each surface based on the elevation, slope, and a curves map analysis and encompass from a few meters to tens of meters, depending on the size of the marker. We checked 80% of the identified surfaces in the field in order to avoid misinterpretation. Some pictures are provided in the Supplement. The dip direction and angle of the surface are computed in a two steps approach. First, we fit a plan using extracted points from the DEM inside the delimited surface. Second, based on this plan we remove the DEM points with residuals 3 times larger than the standard error and compute more accurate plan parameters (second fitting). This outlier suppression removes any inaccurate DEM points and corrects for inaccurate surface delimitation (e.g., integration of a part of the edge of a strath terrace or diffusion processes marks).

Because no obvious initially horizontal markers are known, we propose to correct the marker current slope by the initial one that quantifies the tilt since the marker displacement. To do so we follow the method used by Champagnac et

al. (2008) for the ForeAlps (Western Alps, France). We identify the drain related to the marker formation and compute its current local slope and direction. This method assumes that landscapes are at the equilibrium state and that the river slope has remained constant since the marker formation. This assumption seems reasonable given the major river profiles and because most of the markers used are far from the watershed high-altitude areas, precluding a regressive erosional signal. Finally, we removed the local river plan from the DEM extracted surface.

## 5.2 Morphometrical results

Following this methodology, we obtained 61 surfaces (e.g., strath terraces). We then applied three quality criteria to ensure the robustness of our results, which are as follows: (1) the minimal surface considered is 2500 m², based on a comparison between the 5 m resolution DEM and a real-time kinematic (RTK) GPS survey over three strath terraces (Hérault River); (2) final plans with dip angles larger than 2° are removed; and (3) the residuals for each geomorphological marker must be randomly distributed without marker edge signal or clear secondary structuration. Only 38 markers meet those three quality criteria.

If the identified and corrected markers have indeed registered a differential uplift between the north and the south, we expected the following signals:

- The dipping direction of the tilted markers should be parallel to the main gradient of the topography, i.e., between 150 and 180° E for our studied region. This expectation is the most important one, with regard to uncertainties in the uplift rate and lithospheric elastic parameters.

- There should be a latitudinal tilting trend, i.e., an increase in the tilt angle along the topography gradient. Indeed, null or small tilts are expected near the shoreline and within the maximum uplift area of the Cévennes and Massif Central, while the maximum tilt is expected at a mid-distance between these two regions, i.e., about 50 km inland from the shoreline.

- There should be a positive altitudinal tilting trend (an increase in dip angle with altitude). This trend would be representative of the accumulation of finite tilt. However, it supposes a linear relationship between the altitude and the age of the marker formation. If, at first order, this straightforward hypothesis seems reasonable for river-controlled markers (e.g., strath terraces), other surfaces are hardly expected to follow such an easy relationship.

Among the three expected signal, southward dipping is robustly recorded with a mean tilt angle of 0.60 ± 0.40° with an azimuth of 128° N ± 36° E (Fig. 10). Latitudinal trend and altitudinal trend are less robustly reached, but that is not
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Tilting map of geomorphological benchmark (blue areas). Base map is 30 m resolution DEM with slope shadow. Arrows are oriented according to the marker downward dip. The arrow size is set according to the corrected tilting angle (the bigger the arrow, the more the tilting). Yellow and brown arrows are for robust and less robust surfaces, respectively. Several arrows are hidden because of their small size and too high of proximity to bigger ones. Enumerated site (1) is the Rieutord canyon, (2) is the Leicasse cave system, (3) is the Garrel cave system, and (4) is the Lodève basin with dated basaltic flows. See Fig. 1 for geographical coordinates.</figcaption></figure>
surprising because of the strong susceptibility to local phenomenon or, even so, a lack of robust age constraint.

## 6 Discussion

Both geomorphological and geochronological evidence suggest a Pliocene–Quaternary uplift of the Cévennes area. The origin of such an uplift could be associated with several processes: erosion-induced isostatic rebound, dynamic topography due to mantle convection, thermal isostasy, residual flexural response due to the Gulf of Lion formation, etc. For the Alps and Pyrenees mountains, isostatic adjustment due to erosion and glacial unloading has been recently quantified (Champagnac et al., 2007; Vernant et al., 2013; Genti, 2015; Chéry et al., 2016). Because the erosion rates measured in the Cévennes are similar to those of the eastern Pyrenees (Calvet et al., 2015; Sartégou et al., 2018a), we investigate by numerical modeling how an erosion-induced isostatic rebound could impact the southern Massif Central morphology and deformation.

We define a representative cross section parallel to the main topographic gradient (i.e., NNW–SSE) and close to the field investigation areas (Fig. 11). We study the lithospheric elastic response to erosion with the 2-D finite element model ADELI (Hassani and Chery, 1996; Chéry et al., 2016). The model is composed of a plate accounting for the elasticity of both the crust and uppermost mantle. Although the lithosphere rigidity of the European plate in southern Massif Central is not precisely known, vertical gradient temperatures provided by borehole measurements are consistent with heat flow values ranging from 60 to 70 mW m⁻² (Lucazeau and Vasseur, 1989). Therefore, we investigate plate thickness ranging from 10 to 50 km as done by Stewart and Watts (1997) for studying the vertical motion of the alpine forelands.

We choose values for Young's and Poisson parameters of 10¹¹ Pa and 0.25, respectively, both commonly used values for lithospheric modeling (e.g., Kooi and Cloething, 1992; Champagnac et al., 2007; Chéry et al., 2001). This leads to long-term rigidity of the lithosphere model ranging from 10²¹ to 10²⁵ N m. Since the effect of mantle viscosity on elastic rebound is assumed to be negligible at the timescale of our models (1 to 2 Myr), we neglect the viscoelastic behavior of the mantle. Therefore, the base of the model is supported by an hydrostatic pressure boundary condition balancing the weight of the lithosphere (Fig. 11). Horizontal displacements on vertical sides are set to zero since geodetic measurements show no significant displacements (Nocquet and Calais, 2003; Nguyen et al., 2016). The main parameters controlling our model are the erosion (or sedimentation) triggering isostatic rebound and the elastic thickness.

The erosion profile (Fig. 11) is based on topography, our newly proposed incision rate, and other studies (Olivetti et al., 2016, for onshore denudation and Lofi et al., 2003; Ler-

oux et al., 2014, for offshore sedimentation). This profile is a simplification of the one that can be expected, from Olivetti et al. (2016), and does not aim to precisely match the published data because of the following: (i) the explored time span (~1 Myr) is not covered by thermochronological data (>10 Myr) or cosmogenic denudation rate (10s–100s kyr); (ii) we base our erosion rate on being linked with local (10s km²) slopes that are higher near the drainage divide. We, by this aim, can invoke any kind of erosion processes (e.g., landslides); and (iii) the model assumes a cylindrical structure, and consequently, high-frequency lateral variations in term or actual denudation rate or proxy (slope, elevation, etc.) must be averaged. Concerning this erosion profile, parametric study (highest erosion rate ranging from 1 to 1000 m Ma⁻¹) gives no difference in the interpretation, and for a few percent change in the used erosion rate, the model result (uplift rate) changes by a few percent.

The flexural rigidity controls the intensity and wavelength of the flexural response and ranges from 10²¹ to 10²⁵ N m. It can be expressed as a variation in elastic thickness (T_e) ranging from 4.4 to 96 km (Fig. 12). We also test a possible T_e variation between inland and offshore areas. For the following discussion, we use an elastic thickness of 15 km corresponding to a value of D (flexural rigidity) of 3.75 × 10²³ N m. In this case, the inland and offshore parts are largely decoupled and the large sedimentation rate in the Gulf of Lion does not induce a flexural response on the Cévennes and Grands Causses areas.

With a maximum erosion rate of 80 m Ma⁻¹ (Fig. 11), the models display uplift rates of 50 m Ma⁻¹ over more than 100 km. As previously explained, the finite incision is permitted by an equal amount of uplift considering that the incision is not due to regressive erosion.

Every model shows a general uplift. However, the uplift amplitudes are smaller than the expected ones. To obtain the same uplift rate as the incision rates, the applied erosion rate over the model must be increased. However, we assume that the landscape is at equilibrium, so, if the erosion rate is increased, it will be higher than the incision rate leading to the decay of relief over the area. No evidence of such evolution is found over the region, and if further studies need to be done to quantify the actual erosion rate, we mostly think that a second process is acting, inducing the rest of the uplift that cannot be obtained by the erosion-induced isostatic adjustment. Finally, models predict a seaward tilt of the surface at the regional scale (Fig. 13), in agreement with the observed tilting of morphological markers.

We assume that the sediments collected in the karst were deposited per descensum; i.e., we do not know if the galleries existed a long time before or were formed just before the emplacement of the sediments, but the more elevated the sediments are, the older their deposit is. If there is no evidence of an important aggradation episode leading to more a complex evolution as proposed for the Ardèche canyon (Moccochain et al., 2007; Tassy et al., 2013), we point out that
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Tilting and azimuth distribution. Panel (a) is density distribution for surface maximum tilting in degree. Panel (b) is azimuth of maximum dipping relative to the north. For each histogram, red and grey populations are for robust and primary detected markers.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** (a) Schematic topographic profile. The red box delimits the area shown Figs. 1 and 9. (b) Surface processes profile; negative values are for erosion, and positive values are for sedimentation. (c) Model setup with two compartments (one for the Cévennes area and the second one for the Gulf of Lion). The base of the model is compensated in pressure, and the right and left limits are fixed at zero horizontal velocity and free vertical velocity. $T_c$ is the equivalent elastic thickness (in kilometers), $E$ (in pascals) and $v$ are the Young's modulus and the Poisson coefficient, respectively, whose values are independent in each compartment.</figcaption></figure>

small aggradation or null erosion period could, however, be possible. Some processes could explain such a relative stability, e.g., variation in erosion (due to climatic fluctuation) or impact of eustatic variations (in river profile, flexural response, etc.). Such transient variations have been shown for the Alps (Saillard et al., 2014; Rolland et al., 2017) and are

proposed as being related to climato-eustatic variations and therefore should last 10 to 100 kyr at most.

Based on our sampling resolution, we cannot evidence such transient periods, and we must use an average base level lowering rate in the karst, which we correlate to the incision of the main rivers. The TCN-based incision rate de-
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Modeled uplift according to different $T_e$. Most plausible $T_e$ are between 10 and 30 km.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Modeling result for $T_e = 15$ km. Erosion–sedimentation rate profile is the same as in Fig. 6. Velocity field is shown using arrow for orientation; velocity magnitudes are quantified by the font color code. Black values on top are distance relative to the seashore (positive value landward and negative values seaward). Red line represents the southward, modeled tilting due to differential uplift.</figcaption></figure>

rived from the Rieutord samples ($83 \pm 35$ m Ma⁻¹) is consistent with the one derived from the Garrel samples (U–Th ages: 85.83 m Ma⁻¹, according to the sole U–Th exploitable result; Camus, 2003) and from the Garrel–Leicasse combination (paleomagnetic approach: $84 \pm 21/ - 12$ m Ma⁻¹).

This mean incision rate of ca. $85$ m Ma⁻¹ lasting at least 4 Ma, highlights the importance of the Pliocene–Quaternary period in the Cévennes and Grands Causses morphogenesis. Furthermore, the 300 to 400 m of incision precludes a relative base level controlled by a sea-level drop. Indeed, documented sea level variations are less than 100 m (Haq, 1988; Miller et al., 2005). Furthermore, the Hérault River does not show any significant knickpoints or evidence of unsteadiness in its profile as would be expected if the incision was due to eustatic variations. Therefore, we propose that the incision rate of ~$85$ m Ma⁻¹ is due to a Pliocene–Quaternary uplift of the Cévennes and Grands Causses region.

Other river–valley processes could lead to a local apparent high incision rate, for instance a major landslide or alluvial fan (Ouimet et al., 2008). This hypothesis of an epigenetic formation of the Rieutord is irrelevant because (i) none of the possible causes had been found in the Rieutord canyon and (ii) the consistency of the TCN-based incision rate and the paleomagnetic-based incision rate for two other cave systems. Indeed, the use of two independent approaches and three locations is a good argument in favor of the robustness of our proposed mean $85$ m Ma⁻¹ inci-

sion rate. Yet, using more data, particularly burial dating co-localized with clay samples and adding sampling sites would give a stronger statistical validation. In the Lodève basin (Point 4, Fig. 1), inverted reliefs allow another independent way to quantify minimal incision rate. The K/Ar and paleomagnetic dated basaltic flows spanning from 1 to 2 Myr old that were deposited at the bottom of the former valley (Dautria et al., 2010) are now located at ca. 150 m above the current riverbed, leading to an average incision rate of $77 \pm 10$ m Ma⁻¹, in agreement with karst-inferred incision rates.

Furthermore, preliminary results from canyons on the other side of the Grands Causses (Tarn and Jonte) based on in situ terrestrial cosmogenic dating suggest similar incision rates (Sartégou et al., 2018b) and confirm a regional base level lowering of the Cévennes and Grands Causses region during the Pliocene–Quaternary. This is consistent with the similarities of landscapes and lithologies observed both on the Atlantic and Mediterranean watersheds (e.g., Tarn river).

Once the regional pattern of the Pliocene–Quaternary incision established for the Cévennes–Grands Causses area, the next question is how this river downcutting is related to the regional uplift? First order equilibrium shape and absence of major knick points in the main river profiles preclude the hypothesis of regressive erosion. Hence, back to the three conceptual models presented in part 1 (Fig. 2); we can discard, at first order, the models A (old uplift–old incision) and
B (old uplift-recent incision) because the obtained incision
rate shows recent incision and surface tilting tends to prove
a current uplift. Therefore, the incision rate has to be bal-
anced to the first order by the uplift rate. We add that eustatic
variations are too low in magnitude (100-120 m) and cannot
explain such total incision (up to 400 m). Furthermore, no
obvious evidence of active tectonics is reported for the area,
raising the questions about the processes responsible for this
regional uplift. Very few denudation rates are reported for
our study area (Schaller et al., 2001; Molliex et al., 2016;
Olivetti et al., 2016), and converting canyon incision rates
into denudation and erosion rates is not straightforward, es-
pecially given the large karst developed in the area. Using
a first order erosion-sedimentation profile, that follows the
main topography gradient direction we have modeled the
erosion-induced isostatic rebound. If this process could cre-
ate between half and two-thirds of the Pliocene-Quaternary
uplift, a previously existing topography is needed to trigger
erosion, so it cannot explain the onset of the canyon-carving
or the full uplift rates. Other, processes have to be explored
such as dynamic topography or thermal anomaly beneath the
Massif Central; the magmatism responsible for the impor-
tant increase in volcanic activity since ~6 Myr (Michon and
Merle, 2001; Nehlig et al., 2003) could play a major role,
notably in the initiation of Pliocene-Quaternary uplift. Fur-
ther studies should aim to address the problem of uplift onset,
giving more clues concerning the stable continental area, but
owing to the data we presently have, discussing such an onset
is out of the scope of this paper.

7 Conclusion

Main results of this study are the following three points:

1. Mean incision rate of the Cévennes area is 83 +17/ - 5 m Ma⁻¹ during the last 4 Ma.

2. This incision is due to regional uplift with higher vertical velocities northward.

3. This uplift is partly due (1/2 to 2/3) to isostatic adjust-
ment induced by erosion.

Furthermore, our study highlights the importance of mul-
tidisciplinary approach especially in the study of low-
deformation rate areas.

To the contrary of previous studies that focused on one
cave, we have shown that combining karst burial ages and
paleomagnetic analysis of clay deposits in several caves over
a large elevation range can bring good constraints on incision
rates. This multi-cave system approach diminishes the intrin-
sic limits of the two single methods: low sampling density
(and analysis cost) for the TCN ages and a difficulty setting
the position of paleomagnetic results. Our estimated paleo-
base-level ages are Pliocene-Quaternary (ca. last 4 Ma) and
allow the derivation of a mean incision rate of 83 +17/ -

5 m Ma⁻¹ for the Cévennes area. The landscape and espe-
cially the river profiles suggest a first-order equilibrium, al-
lowing the consideration of the incision rate as an uplift rate.

We have shown using a geomorphological analysis that at
least south of the Cévennes, several surfaces are tilted to-
ward the SSE. This kind of study had been performed before
on large structures (Champagnac et al., 2007) or endokarstic
markers (Granger and Stock, 2004), but it is the first time that
it is performed at such a scale with small markers. Numerical
modeling yields the same pattern of SSE dipping, allowing
more confidence in the geomorphometric results.

Our multidisciplinary approach brings the first absolute
dating of the Cévennes landscapes and suggests that the
present-day morphology is partly inherited from the Plio-
Quaternary erosion-induced isostatic rebound.

We propose that related erosional isostatic adjustment is
of major importance for the understanding of the southern
French Massif Central landscape evolution and explains a
large part of the uplift.

At a larger scale, we assume that the main conclusion of
our study can be extrapolated to the majority of the intraplate
orogens. That is to say, once the forces responsible for the
initial uplift (e.g., plate tectonics or dynamic topography)
fade out, the uplift continues thanks to erosion-induced iso-
static adjustment.

An analysis at the scale of the Massif Central is now
needed before nailing down our interpretations of the Mas-
sif Central dynamics.

*Code and data availability.* Surface analysis was performed using QGIS version 2.18, MATLAB® code and IGN DEM (RGE Alti®) 5 m. Modeling was performed using ADELI code (Hassani and Chery, 1996; Chéry et al., 2016). Data for TCN and paleomagnetic analysis are provided in the paper itself or in the Supplement. Additional information for geologic background are available at http://infoterre.brgm.fr/ (last access: 21 February 2020) (French Geological Survey data visualizer) and in Karst3D Team (2019, https://doi.org/10.15148/940c2882-49f1-49db-a97e-12303cace752).

*Supplement.* The supplement related to this article is available on-line at: https://doi.org/10.5194/se-11-241-2020-supplement.

*Author contributions.* OM, PV, and GC did the sampling. GC and DF performed the TCN analysis. PC and OM did the magnetic measurements and interpretations. OM did the surface identification and analysis. OM, PV, and JC performed the numerical model. OM, PV, JFR, GC, PC, JC and DF interpreted and wrote the paper.

***Competing interests.*** The authors declare that they have no conflict of interest.
*Acknowledgements.* We are grateful to ANSTO for providing facilities for chemical extraction for the TCN analysis. We thank the reviewers for useful remarks and comments that we think help to increase the level of the paper.

*Review statement.* This paper was edited by Federico Rossetti and reviewed by two anonymous referees.

## References

Arthaud, F. and Laurent, P.: Contraintes, déformations et déplacements dans l'avant-pays pyrénéen du Languedoc méditerranéen, Godin. Acta, 8, 142-157, 1995.

Audra, P., Camus, H., and Rochette, P.: Le karst des plateaux de la moyenne vallée de l'Ardèche: datation par paléomagnétisme des phases d'évolution plio-quaternaires (aven de la Combe Rajeau), Bull. Soc. Géol. France, 172, 121-129, 2001.

Balco, G., Stone, J. O., Lifton, N. A., and Dunai, T. J.: A complete and easily accessible means of calculating surface exposure ages or erosion rates from Be-10 and Al-26 measurements, Quat. Geochronol., 3, 174-195, 2008.

Barbarand, J., Lucazeau, F., Pagel, M., and Séranne, M.: Burial and exhumation history of the south-eastern Massif Central (France) constrained by en apatite fission-track thermochronology, Tectonophysics, 335, 275-290, 2001.

Barruol, G. and Granet, M.: A Tertiary asthenospheric flow beneath the southern French Massif Central indicated by upper mantle seismic anisotropy and related to the west Mediterranean extension, Earth Planet. Sci. Lett., 202, 31-47, 2002.

Brichau, S., Respaut, J. P., and Monié, P.: New age constraints on emplacement of the Cévenol granitoids, South French Massif Central, Int. J. Earth Sci. 97, 725-738, https://doi.org/10.1007/s00531-007-0187-x, 2007.

Bruxelles, L.: Dépôts et altérites des plateaux du Larzac central: causses de l'Hospitalet et de Campestre (Aveyron, Gard, Hérault) Evolution morphogénétique, conséquences géologiques et implications pour l'aménagement, Thèse, spécialité: Miliéux physiques méditerranéens, Université d'Aix-Marseille I, Université de Provence, UFR Sciences géographiques et de l'aménagement, 2001.

Calais, E., Freed, A. M., Van Arsdale, R., and Stein, S.: Triggering of New Madrid seismicity by late-Pleistocene erosion, Nature, 466, 608-611, https://doi.org/10.1038/nature09258, 2010.

Calais, E., Camelbeeck, T., Stein, S., Liu, M., and Craig, T. J.: A new paradigm for large earthquakes in stable continental plate interiors, Geophys. Res. Lett., 43, 10621-10637, https://doi.org/10.1002/2016GL070815, 2016.

Calvet, M., Gunnell, Y., Braucher, R., Hez, G., Bourlès, D., Guillou, V., Delmas, M., and ASTER team: Cave levels as proxies for measuring post-orogenic uplift: Evidence from cosmogenic dating of alluvium-filled caves in the French Pyrenees, Geomorphology, 246, 617-633, https://doi.org/10.1016/j geomorph.2015.07.013, 2015.

Camus, H.: Vallée et réseaux karstiques de la bordure carbonatée sud-cévenole. Relation avec la surrection, le volcanisme et les paléoclimats, Thèse de doctorat, Université Bordeaux, Bordeaux, 3, 692 pp., 2003.

Champagnac, J. D., Molnar, P., Anderson, R. S., Sue, C., and Delacou, B.: Quaternary erosion-induced isostatic rebound in the western Alps, Geology, 35, 195-198, https://doi.org/10.1130/G23053A.1, 2007.

Champagnac, J.-D., van der Beek, P., Diraison, G., and Dauphin, S.: Flexural isostatic response of the Alps to increased Quaternary erosion recorded by foreland basin remnants, SE France, Terra Nova, 20, 213-220, https://doi.org/10.1111/j.1365-3121.2008.00809.x, 2008.

Chéry, J., Zoback, M. D., and Hassani, R.: An integrated mechanical model of the San Andreas Fault in central and northern California, J. Geophys. Res., 106, 22051-22066, 2001.

Chéry, J., Genti, M., and Vernant, P.: Ice cap melting and low-viscosity crustal root explain the narrow geodetic uplift of the Western Alps, Geophys. Res. Lett., 43, 1-8, 2016.

Child, D. P., Elliott, G., Mifsud, C., Smith, A. M., and Fink, D.: Sample processing for earth science studies at ANTARES, Nuclear Instruments and Methods in Physics Research Section B Beam Interactions with Materials and Atoms, 172, 856-860, https://doi.org/10.1016/S0168-583X(00)00198-1, 2000.

Corbel, J.: Les phénomènes karstiques dans les Grands Causses, in: Revue de géographie de Lyon, 29, 287-315, https://doi.org/10.3406/geoca.1954.1990, 1954.

Dautria, J. M., Liotard, J. M., Bosch, D., and Alard, O.: 160 Ma of sporadic basaltic activity on the Languedoc volcanic line (Southern France): A peculiar case of lithosphere-asthenosphere interplay, Lithos, 120, 202-222, https://doi.org/10.1016/j.lithos.2010.04.009, 2010.

Genti, M.: Impact des processus de surface sur la déformation actuelle des Pyrénées et des Alpes, Géophysique, Thèse, Université de Montpellier, Montpellier, 2015 (in French).

Granet, M., Wilson, M., and Achauer, U.: Imaging a mantle plume beneath the French Massif Central, Earth Planet. Sc. Lett., 136, 281-296, 1995.

Granger, D. E. and Muzikar, P. F.: Dating sediment burial with in situ-produced cosmogenic nuclides: theory, techniques, and limitations, Earth Planet. Sc. Lett., 188, 269-281, 2001.

Granger, D. E. and Stock, G. M.: Using cave deposits as geologic tiltmeters: Application to postglacial rebound of the Sierra Nevada, California, Geophys. Res. Lett., 31, L22501, https://doi.org/10.1029/2004GL021403, 2004.

Granger, D. E., Kirchner, J. W., and Finkel, R. C.: Quaternary downcutting rate of the New River, Virginia, measured from differential decay of cosmogenic $^{26}$Al and $^{10}$Be in cave-deposited alluvium, Geology, 25, 107-110, 1997.

Granger, D. E., Fabel, D., and Palmer, A. N.: Pliocene-Pleistocene incision of the Green River, Kentucky determined from radioactive decay of cosmogenic $^{26}$Al and $^{10}$Be in Mammoth Cave sediments, GSA Bulletin, 113, 825-836, 2001.

Granger, D. E., Gibbon, R. J., Kuman, K., Clarke, R. J., Bruxelles, L., and Caffee, M. W.: New cosmogenic burial ages for Sterkfontein Member 2 Australopithecus and Member 5 Oldowan, Nature Lett., 522, 85, https://doi.org/10.1038/nature14268, 2015.

Hajna, Z. N., Mihevc, A., Pruner, P., and Bosák, P.: Palaeomagnetic research on karst sediments in Slovenia, Int. J. Speleol., 39, 47-60, 2010.

Haq, B. U., Herdenbol, J., and Vail, P. R.: Mesozoic and cenozoic chronostratigraphy and cycles of sea-level change, Society
Economic Paleontologists Mineralogists Special Publication, 42, 71–108, Tulsa, Oklahoma, 1988.

Harmand, D., Adamson, K., Rixhon, G., Jaillet, S., Losson, B., Devos, A., Hez, G., Calvet, M., and Audra, P.: Relationships between fluvial evolution and karstification related to climatic, tectonic and eustatic forcing in temperate regions, Quat. Sci. Rev., 2017, 1–19, https://doi.org/10.1016/j.quascirev.2017.02.016, 2017.

Hassani, R. and Chery, J.: Anaelasticity explains topography associated with Basin and Range normal faulting, Geology, 24, 1095, https://doi.org/10.1130/0091-7613(1996)024<1095:AETAWB>2.3.CO;2, 1996.

Hill, C. A.: Sedimentology and Paleomagnetism of sediments, Kartchner caverns, Arizona, J. Cave Karst Studies, 61, 79–83, 1999.

Husson, E.: Intéraction géodynamique/karstification et modélisation 3D des massifs carbonatés: Implication sur la distribution prévisionnelle de la karstification. Exemple des paléokarsts créticés à néogènes du Languedoc montpelliérain, Sciences de la Terre, Université Montpellier 2 – Sciences et techniques du Languedoc, 236 pp., 2014.

Karst3D Team: KARST3D, OSU OREME (Collection), https://doi.org/10.15148/940c2882-49f1-49db-a97e-12303cace752, 2019.

Kooi, H., Cloetingh, S., and Burrus, J.: Lithospheric Necking and Regional Isostasy at Extensional Basins 1. Subsidence and Gravity Modeling With an Application to the Gulf of Lions Margin (SE France), J. Geophys. Res., 97, 17553–17571, 1992.

Leroux, E., Rabineau, M., Aslanian, D., Granjeon, D., Droz, L., and Gorini, C.: Stratigraphic simulations of the shelf of the Gulf of Lions: testing subsidence rates and sea-level curves during the Pliocene and Quaternary, Terra Nova, 26, 230–238, https://doi.org/10.1111/ter.12091, 2014.

Lofi, J., Rabineau, M., Gorini, C., Berne, S., Clauzon, G., De Clarens, P., Dos Reis, A. T., Mountain, G. S., Ryan, W. B. F., Steckler, M. S., and Fouchet, C.: Plio-Quaternary prograding clinoform wedges of the western Gulf of Lion continental margin (NW Mediterranean) after the Messinian Salinity Crisis, Mar. Geol., 198, 289–317, https://doi.org/10.1016/S0025-3227(03)00120-8, 2003.

Lucazeau, F. and Vasseur, G.: Heat flow density data from France and surrounding margins, Tectonophysics, 164, 251–258, 1989.

Manchuel, K., traversa, P., Baumont, D., Cara, M., Nayman, E., and Durouchoux, C.: The French seismic CAT-alogue (FCAT-17), Bull. Earthquake Eng., 16, 2227–2251, https://doi.org/10.1007/s10518-017-0236-1, 2018.

Miallier, D., Michon, L., Évin, J., Pilleyre, T., Sanzelle, S., and Vernet, G.: Volcans de la chaîne des Puys (Massif central, France): point sur la chronologie Vasset-Kilian-Pariou-Chopine, Comptes Rendus Geoscience, 336, https://doi.org/10.1016/j.crte.2004.08.002, 2004.

Michon, L. and Merle, O.: The evolution of the Massif Central rift: Spatio-temporal distribution of the volcanism, Bulletin de la Society Geologique de France, 172, 201–211, https://doi.org/10.2113/172.2.201, 2001.

Miller, K. G., Kominz, M. A., Browning, J. V., Wright, J. D., Mountain, G. S., Katz, M. E., Sugarman, P. J., Cramer, B. S., Christie-Blick, N., and Pekar, S. F.: The Phanerozoic record of global sea-level change, Science, 310, 1293–1298, https://doi.org/10.1126/science.1116412, 2005.

Mocochain, L.: Les manifestations geodynamiques – Externes et internes – de la crise de salinité messinienne sur une plate-forme carbonatée peri-méditerranéenne: le karst de la basse ardèche (moyenne vallée du rhône, France), Thèse de doctorat, Université Aix-Marseille I – Université de Provence U.F.R des Sciences géographiques et de l'aménagement Centre Européen de Recherches et d'Enseignement en Géosciences de l'Environnement, 196 pp., 2007.

Molliex, S., Rabineau, M., Leroux, E., Bourlès, D. L., Authemayou, C., Aslanian, D., Chauvet, F., Civet, F., and Jouët, G.: Multi-approach quantification of denudation rates in the Gulf of Lion source-to-sink system (SE-France), Earth Planet. Sc. Lett., 444, 101–115, https://doi.org/10.1016/j.epsl.2016.03.043, 2016.

Nehlig, P., Boivin, P., de Goër, A., Mergoil, J., Prouteau, G., Strac, G., and Thiéblemont, D.: Les volcans du Massif central, Revue BRGM: Géologues, Numéro Spécial: Massif central, Orleans (France), 2003.

Nguyen, H. N., Vernant, P., Mazzotti, S., Khazaradze, G., and Asensio, E.: 3-D GPS velocity field and its implications on the present-day post-orogenic deformation of the Western Alps and Pyrenees, Solid Earth, 7, 1349–1363, https://doi.org/10.5194/se-7-1349-2016, 2016.

Nishiizumi, K., Imamura, M., Caffee, M. W., Southon, J. R., Finkel, R. C., and McAninch, J.: Absolute calibration of 10Be AMS standards, Nucl. Instrum. Meth. B, 258, 403–413, 2007.

Nocquet, J.-M. and Calais, E.: Crustal velocity field of western Europe from permanent GPS array solutions, 1996–2001, Geophys. J. Int., 154, 72–88, https://doi.org/10.1046/j.1365-246X.2003.01935.x, 2003.

Nocquet, J.-M., Sue, C., Walpersdorf, A., Tran, T., Lenôtre, N., Vernant, P., Cushing, M., Jouanne, F., Masson, F., Baize, S., Chéry, J., and Van der Beek, P. A.: Present-day uplift of the western Alps, Sci. Rep.-UK, 6, 28404, https://doi.org/10.1038/srep28404, 2016.

Olivetti, V., Godard, V., Bellier, O., and ASTER team: Cenozoic rejuvenation events of Massif Central topography (France): Insights from cosmogenic denudation rates and river profiles, Earth Planet. Sc. Lett., 444, 179–191, https://doi.org/10.1016/j.epsl.2016.03.049 0012-821X, 2016.

Ouimet, W. B., Whipple, K. X., Crosby, B. T., Johnson, J. P., Schildgen, T. F.: Epigenetic gorges in fluvial landscapes, Earth Surf. Process. Landf., 33, 1993–2009, https://doi.org/10.1002/esp.1650, 2008.

Rolland, Y., Petit, C., Saillard, M., Braucher, R., Bourlès, D., Darnault, R., Cassol, D., and ASTER Team: Inner gorges incision history: A proxy for deglaciation? Insights from Cosmic Ray Exposure dating (¹⁰Be and ³⁶Cl) of river-polished surfaces (Tinée River, SW Alps, France), Earth Planet. Sc. Lett., 457, 271–281, https://doi.org/10.1016/j.epsl.2016.10.007, 2017.

Rovey II, C. W., Balco, G., Forir, M., and Kean, W. F.: Stratigraphy, paleomagnetism, and cosmogenic-isotope burial ages of fossil-bearing strata within Riverbluff Cave, Greene County, Missouri, Quaternary Res., 2017, 1–13, https://doi.org/10.1017/qua.2017.14, 2017.

Saillard, M., Petit, C., Rolland, Y., Braucher, R., Bourlès, D. L., Zerathe, S., Revel, M., and Jourdon, A.: Late Quaternary incision rates in the Vésubie catchment area (Southern French Alps) from record of global sea-level change, Science, 310, 1293–1298, https://doi.org/10.1126/science.1116412, 2005.
in situ-produced ³⁶Cl cosmogenic nuclide dating: Tectonic and climatic implications, J. Geophys. Res.-Earth Surf., 119, 1121–1135, https://doi.org/10.1002/2013JF002985, 2014.

Sanchis, E. and Séranne, M.: Structural style and tectonic evolution of a polyphase extensional basin of the Gulf of Lion passive margin: the Tertiary Alès basin, southern France, Tectonophysics, 322, 219–242, https://doi.org/10.1016/S0040-1951(00)00097-4, 2000.

Sartégou, A.: Évolution morphogénique des Pyrénées orientales: apports des datations de systèmes karstiques étagés par les nucléides cosmogéniques et la RPE, Géomorphologie. Thèse de l'Université de Perpignan, Perpignan, France, 2017 (in French).

Sartégou, A., Bourlès, D. L., Blard, P.-H., Braucher, R., Tibari, B., Zimmermann, L., Leanni, L., Aumaitre, G., and Keddadouche, K.: Deciphering landscape evolution with karstic networks – A Pyrenean case study, Quat. Geochronol., 43, 12–29, https://doi.org/10.1016/j.quageo.2017.09.005, 2018a.

Sartégou, A., Mialon, A., Thomas, S., Giordani, A., Lacour, Q., Jacquet, A., André, D., Calmels, L., Bourlès, D. L., Bruxelles, L., Braucher, R., Leanni, L., and ASTER team: When TCN meet high school students: deciphering western Cévennes landscape evolution (Lozère, France) sin g TCN on karstic networks, Poster 4th Nordic Workshop on Cosmogenic Nuclides, 4–6 June 2018, Geiranger, Norway, 2018b.

Schaller, M., von Blanckenburg, F., Hovius, N., and Kubik, P. W.: Large-scale erosion rates from in situ-produced cosmogenic nuclides in European river sediments, Earth Planet. Sc. Lett., 188, 441–458, 2001.

Séranne, M., Benedicto, A., Labaum, P., Truffert, C., and Pascal, G.: Structural style and evolution of the Gulf of Lion Oligo-Miocene rifting: role of the Pyrenean orogeny, Mar. Petrol. Geol., 12, 809–820, 1995.

Séranne, M., Camus, H., Lucazeau, F., Barbarand, J., and Quinif, Y.: Surrection et érosion polyphasées de la Bordure cévenole. Un exemple de morphogenèse lente, Bull. Soc. Géol. France, 173, 97–112, 2002.

Sibuet, J.-C., Srivastava, S. P., and Spakman, W.: Pyrenean orogeny and plate kinematics, J. Geophys. Res.-Sol. Ea., 109, B08104, https://doi.org/10.1029/2003JB002514, 2004.

Spassov, S. and Valet, J.-P.: Detrial magnetisations from redeposition experiments of different natural sediments, Earth Planet. Sc. Lett., 351–352, 147–157, https://doi.org/10.1016/j.epsl.2012.07.016, 2012.

Stewart, J. and Watts, A. B.: Gravity anomalies and spatial variation of flexural rigidity at mountain ranges, J. Geophys. Res., 102, 5327–5352, https://doi.org/10.1029/96JB03664, 1997.

Stock, G. M., Granger, D. E., Sasowsky, I. D., Anderson, R. S., and Finkel, R. C.: Coomparison of U-Th, paleomagnetism, and cosmogenic burial methods for dating caves: Implications for landscape evolution studies, Earth Planet. Sc. Lett., 236, 388–403, https://doi.org/10.1016/j.epsl.2005.04.024, 2005.

Tarayoun, A., Mazzotti, S., and Gueydan, F.: Quantitative impact of structural inheritance on present-day deformation and seismicity concentration in intraplate deformation zones, Earth Planet. Sc. Lett., 518, 160–171, https://doi.org/10.1016/j.epsl.2019.04.043, 2017.

Tassy, A., Mocochain, L., Bellier, O., Braucher, R., Gattacceca, J., and Bourlès, D.: Coupling cosmogenic dating and magnetostatigraphy to constrain the chronological evolution of peri-Mediterranean karsts during the Messinian an the Pliocene: Example of Ardèche Valley, Southern France, Geomorphology, 189, 81–92, https://doi.org/10.1016/j geomorph.2013.01.019, 2013.

Tauxe, L., Steindorf, J. L., and Harris, A.: Depositional remanent magnetisation: Toward an improved theatrical and experimental foundation, Earth Planet. Sc. Lett., 244, 515–529, https://doi.org/10.1016/J.epsl.2006.02.003, 2006.

Tricart, P.: From passive margin to continental collision: A tectonic scenario for the western Alps, Am. J. Sci., 284, 97–120, 1984.

Vernant, P., Hivert, F., Chéry, J., Steer, P., Cattin, R., and Rigo, A.: Erosion-induced isostatic rebound triggers extension in low convergent mountain ranges, Geology, 41, 467–470, https://doi.org/10.1130/G33942.1, 2013.