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

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 (¹⁰Be∕²⁶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 $+\mathrm{17}/-\mathrm{5}$ m Ma⁻¹ 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³–10⁵ 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 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; Barbarand 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 ¹⁰Be and ²⁶Al isotopes. The ¹⁰Be∕²⁶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.

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.

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érrane et al., 2002), we would find surface tilting, but cosmogenic burial dating with ¹⁰Be∕²⁶Al, which cannot discern ages older than ∼5 Ma due to excessive decay of ²⁶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 pre-date 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 radiation 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.

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

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). 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 ¹⁰Be and ²⁶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.

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.

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.

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

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

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⁻¹. 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 of 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 totalizing 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.

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

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

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

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 normal 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±9 m Ma⁻¹.

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 m Ma⁻¹ incision rate range (linear rate, each 1 m Ma⁻¹). Modeled polarities are found using the intersection between sample elevation and incision rate line.

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.

We obtained 10 possible incision rates with the same best correlation factor (Fig. 8) spanning from 43 to 111 m Ma⁻¹ (mean of 87±24 m Ma⁻¹). Taking into account the transitional signal of the LGP site in the Leicasse cave system yields a linear incision rate of 83 $+\mathrm{17}/-\mathrm{5}$ m Ma⁻¹. 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±11 m Ma⁻¹ (Fig. 8). Those two results, based on independent computations, suggest the same first-order incision rate for the last 4 Ma of 84 $+\mathrm{21}/-\mathrm{12}$ m Ma⁻¹. 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, inducing 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.

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.

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 address 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 emplacement. 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 recessive 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 surprising because of the strong susceptibility to local phenomenon or, even so, a lack of robust age constraint.

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.

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.

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

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

Figure 12 Modeled uplift according to different T_e. Most plausible T_e are between 10 and 30 km.

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.

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.

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 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 derived from the Rieutord samples (83±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 $+\mathrm{21}/-\mathrm{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⁻¹ incision 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±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 balanced 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, especially 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 create 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 important 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. Further 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 $+\mathrm{17}/-\mathrm{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 adjustment induced by erosion.

Furthermore, our study highlights the importance of multidisciplinary 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 intrinsic 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 $+\mathrm{17}/-\mathrm{5}$ m Ma⁻¹ for the Cévennes area. The landscape and especially the river profiles suggest a first-order equilibrium, allowing 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 toward 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 isostatic adjustment.

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