Cretaceous and late Cenozoic uplift of a Variscan Massif: The case of the French Massif Central studied through low-temperature thermochronometry

Many mountain belts on Earth are not associated with convergent domains and crustal shortening. These intraplate orogens are characterized by low seismic activity, slow deformation rates, and a typically long-wavelength geomorphological pattern with slowly eroding high-elevated plateaus and erosion mainly focused along their margins (e.g., Stanley et al., 2013; Scotti et al., 2014). The Massif Central in France is a striking example of such a nonconvergent mountain; it is located outside of the European-African plate boundary and was not involved in its Cenozoic convergence. It is at present day a prominent topographic feature reaching up to 1700 m above sea level (asl; Fig. 1) with a low-relief and high-elevation landscape bordered by incised fluvial valleys (up to 1000 m of relief; Figs. 1B and 1C), but the processes responsible for the creation and persistence of this relief have not yet been clearly elucidated.

The French Massif Central is not unique, and other portions of the Variscan belt in Central Europe exhibit a present-day high topography associated with a potential Cenozoic rejuvenation. For these massifs, the common proposed uplift phases are interpreted to have been mainly controlled by the far-field effect of the Europa-Africa convergence, such as a shortening event during the Late Cretaceous–Paleocene and Eocene–Miocene intracontinental rifting (Ziegler, 1994; Kley and Voigt, 2008). Moreover, a Cenozoic uplift event has been recognized all over Central Europe and is usually considered to have been associated with volcanism and high-temperature mantle anomalies (Ziegler and Dèzes, 2007). This event occurred with different timing and amplitude, suggesting a complex regional variability that reflects the complexity of the European geodynamic setting.

Important insights on the Cenozoic geodynamics of the Massif Central can be gained from the reconstruction of the timing, rates, and amount of relief growth. Different scenarios can be proposed for the origin of this present-day topography, spanning from (1) a long-lasting persistence of a Variscan mountain range characterized by slow, erosion-induced isostatic uplift or (2) a recent topographic rejuvenation, which is supported by many geomorphological observations, such as recent river incision, multiple levels of perched surfaces (Séranne et al., 2002; Olivetti et al., 2016), and a complex history of burial and exhumation phases highlighted by thermochronology (Barbarand et al., 2001; Peyaud et al., 2005; Gautheron et al., 2009). In the framework of this rejuvenated topography hypothesis, the ages and mechanisms for a regional uplift are not explicitly defined and could encompass different phases, such as a post-Miocene uplift event induced by the Messinian salinity crisis (Mocochain et al., 2009; Tassy et al., 2013), a Miocene uplift phase triggered by mantle upwelling and associated with extensive intraplate volcanism, or an Oligocene event associated with rifting and opening of the Liguro-Provencal Basin (Fig. 1A; Séranne et al., 2002; Faccenna et al., 2010).

The actual timing and mechanisms for the formation of the topography of the Massif Central remain unclear, but the general current interpretation implies a contribution of mantle upwelling leading in turn to volcanic activity that deeply shaped the massif landscape since the Miocene (i.e., Michon and Merle, 2001). Numerous geophysical campaigns were carried out in the 1990s and highlighted a seismic low-velocity zone 150–300 km deep, high heat flux, and negative Bouguer gravity anomaly (Granet et al., 1995; Sobolev et al., 1997). All these observations are consistent with the interpretation of a shallow mantle and thus potential mantle contribution to the late-stage topographic evolution of the Massif Central.

Although several studies have reported geomorphological evidence for river incision throughout the massif (for a review, see Séranne et al., 2002), an attempt to reconstruct the large-scale and long-term topographic evolution of the massif is still lacking. To address this question, we performed a new thermochronological study based on apatite (U-Th)/He and fission-track analysis and integrated the results with published data (Barbarand et al., 2001; Gautheron et al., 2009). Sampling was performed along the eastern margin of the Massif Central, from the high-elevation surfaces to the bottom of the Rhône River valley, in order to cover the entire margin profile (total elevation transect of 1200 m; Fig. 1D). The eastern margin should have recorded the erosion signal induced by any potential uplift within the entire massif because the long persistence of the Rhône River formed a stable regional base level for erosion processes.

The French Massif Central is mainly composed of rocks associated with the Variscan orogeny (Devonian to late Carboniferous), dominated by granites and high-grade metamorphic rocks with minor mafic and ultramafic rocks and late Paleozoic sedimentary cover (Fig. 2A; Chantraine et al., 2003). Paleozoic rocks are overlain by thin Triassic continental sandstones and by Triassic to Cretaceous mainly carbonate sediments that characterized the sedimentation of the entire northern margin of the Tethys basin. During the Late Cretaceous, the Pyrenean orogeny involved mainly the southern portion of the massif, and it reactivated Paleozoic structures with mainly strike-slip kinematics (Blès et al., 1989).

The European Cenozoic rift system event (i.e., Ziegler 1994), recognized throughout Western and Central Europe, involved mainly the northern and eastern parts of the Massif Central, producing N-S–elongated basins filled by several hundreds of meters, and locally up to 1–2 km, of marine to lacustrine sediments (Limagne, Roanne, Bresse, and Valence basins; Figs. 2A and 2B). Along the southern margin of the massif, the European Cenozoic rifting event overlapped with the opening of the Gulf of Lion, which induced extensional tectonics along reactivated Paleozoic features such as the Nimes and Cévennes faults (Figs. 2A and 2B).

The volcanic history of the Massif Central starts back in the early Cenozoic, and it evolved in three main stages (Michon and Merle, 2001, and references therein): (1) a prerift Paleocene phase, with limited magmatic production, (2) a rift-related phase found in the northern part of the massif only, and (3) the main phase from ca. 27 Ma to the Quaternary. During this last phase, the main magmatic activity occurred from the middle Miocene (ca. 15 Ma) to the Quaternary with a climax between 9 and 6 Ma, when the volume of erupted magma was much larger than during any other magmatic event in this area. The erupted magmas are subalkaline to alkaline intraplate-type lavas (Lustrino and Wilson, 2007).

The Alpine orogeny only marginally affected the Massif Central, and although the Alpine frontal thrusts are currently only a few tens of kilometers to the east, no clear evidence of Alpine deformation is reported within the crystalline basement of the massif (Blés et al., 1989).

The Early Triassic topography of the Massif Central, at the end of the Variscan orogenic cycle, is interpreted to have been characterized by low elevations, limited relief, and marginal plains (Le Griel, 1988) bordered by marine basins that hosted carbonate sedimentation (Curnelle and Dubois, 1986). The Variscan basement should have been close to the base level during the whole Jurassic, without evidence of prominent topography capable of producing large fluvial systems and clastic sediments (Le Griel, 1988). In contrast, an important phase of surface uplift and consequent crustal exhumation occurred during the Cretaceous, based on low-temperature thermochronological data (Barbarand et al., 2001; Peyaud et al., 2005; Gautheron et al., 2009) and stratigraphic constraints (Curnelle and Dubois, 1986). This Cretaceous exhumation was coeval with the opening of the Bay of Biscay and slightly predated the emersion of the “Ithsme Durancien,” an aerially exposed E-W–elongated area that emerged during the Cenomanian and led to the formation of large bauxite deposits. The emersion of the “Ithsme Durancien” was associated with a NE-SW–directed extensional tectonic regime that formed the Vocontian Basin to the north and South Provence Basin to the south (Tavani et al., 2018, and references therein). During the Early Cretaceous, topography of the southern Massif Central should have been characterized by long wavelengths of ∼500 km with a moderate amplitude of 300–500 m (Wyns and Guillocheau, 1999) produced by a far-field response to the Pyrenean orogeny. Cretaceous low topography is also suggested by the limited clastic sedimentation within the subsiding marine basins, such as the Grand Causse, coeval with the exhumation of the Cévennes and Ardèche regions (Barbarand et al., 2001; Gautheron et al., 2009).

During the Paleocene, large portions of the southeastern massif, located south of the Cévennes fault, were emerged, as a consequence of the Pyrenean compression. Farther north, where there is not significant evidence of any Pyrenean activity, the topography was probably limited to a few hundred meters (Séranne et al., 2002, and references therein).

Stratigraphic evidence suggests that during the early Eocene, a large portion of the Massif Central was close to sea level. Marine clays and evaporitic deposits are described in the Velay region, in the central part of the massif (Girod et al., 1979; Rey, 1971; Turland, et al., 1994), while the sedimentation in the Bresse basin recorded marine communication with the Valence basin and the Mediterranean Sea (Sissingh, 1998). In the southern part of the massif, marine deposits are described in the Ales basin (Alabouvette and Cavelier, 1984).

European Cenozoic rifting produced deep grabens, mainly in the northern massif (e.g., Limagne basin), with the deposition of up to 2–3 km of sediments characterized by lacustrine-continental facies with regular marine incursions. Rifting induced medium to high crustal thinning (between 5% and 25%) but was not associated with a large emission of volcanic products (Michon and Merle, 2001). The lack of synchronous thermochronological cooling ages also suggests that there was no significant exhumation event (Fig. 2C; Barbarand et al., 2001; Peyaud et al., 2005; Gautheron et al., 2009).

In the center of the massif (Margeride and Velay regions), geometric relationships between river incision and volcanic flows enabled multiple river incision events to be constrained between 10 Ma and 7 Ma, which were attributed to an uplift phase of around 200–300 m (Michon, 2000, and references therein; Goër de Herve and Etienne, 1991; Defive and Cantagrel, 1998), while in the Causse region (southern massif), the incision of the Tarn River was active already at 13 Ma (Ambert, 1994). In general, the Miocene to present incision history of most of the Massif Central area developed first with a minor phase from 13 to 5 Ma, followed by a major phase of river incision, evaluated at around 400 m, suggesting an important and coeval uplift event (Etienne, 1970; Defive and Cantagrel, 1998; Michon, 2000; Séranne et al., 2002). This second uplift phase has been recently detected by river profile analysis and cosmogenic nuclide–derived denudation rates along the eastern margin of the massif, and it has been consistently evaluated to be of the order of 300–400 m (Olivetti et al., 2016).

Previous thermochronological studies have focused along the south and southeastern portions of the Massif Central (Fig. 2C). The Ardèche and Cévennes areas, where the relief is greater, were investigated using apatite fission-track (AFT; Barbarand et al., 2001) and (U-Th)/He (AHe; Gautheron et al., 2009) dating on the same set of samples. AFT and AHe ages were between ca. 130 Ma and 45 Ma, and thermal modeling suggested an exhumation event with onset at ca. 130 ± 10 Ma (Barbarand et al., 2001). The positive correlation between age and elevation for the AHe and AFT data was interpreted to be a consequence of prolonged exhumation since 130 Ma (apparent exhumation rate between 0.01 and 0.02 km/m.y.), leading to slow cooling through the apatite partial annealing zone (PAZ) and the He partial retention zone (PRZ).

A very similar thermal history has been proposed by Peyaud et al. (2005) using AFT ages (spanning between 147 Ma and 77 Ma) and vitrinite reflectance data for the region of the Rouergue and Montagne Noire (Fig. 2C). In addition to this pre-Cenozoic thermal history, Barbarand et al. (2001) observed a characteristic shape of the mean track length (MTL) versus AFT age relationship in the Cévennes mountains, where the youngest and oldest ages display relatively long mean track lengths, while the intermediate ages are associated with shorter track lengths. Such a relationship (“boomerang” shape) is likely the consequence of a period of thermal annealing, affecting a set of samples to various degrees, followed by a discrete cooling event (Green, 1986; Omar et al., 1989). The oldest ages represent shallower and thus cooler rocks, preserving most of their tracks and previous thermal histories, while intermediate ages (the middle “concave-up” section) represent rocks that resided a long time within the PAZ and that are characterized by more severely shortened tracks. The timing of the cooling event is recorded by “deepest” youngest samples, provided that they have MTLs of the order of 14 µm (i.e., they were at paleotemperatures >∼110–125 °C prior to the last phase of cooling); otherwise, they provide just an upper limit to this timing. In the data set of Barbarand et al. (2001), the younger part of the “boomerang” pattern is less well defined, and the youngest sample displayed an age of 45 ± 2 Ma accompanied by a MTL of 12.9 ± 0.2 µm. Additional data are required to further elucidate the recent exhumation history of the Massif Central and its relationship with geodynamic events and potential topographic rejuvenation mechanisms.

In the context of ancient, slowly eroding orogens, thermochronological methods with low closure temperatures offer the best opportunity to detect any recent topographic rejuvenation. For that purpose, we mainly used apatite (U-Th)/He dating for its low closure temperature (40–120 °C; Gautheron et al., 2009; Flowers et al., 2009; Djimbi et al., 2015), complemented by some AFT data (closure temperature of 110° ± 10 °C; Green et al., 1989) for a subset of samples to expand the time window of investigation. Samples were collected along pseudo-elevation profiles, tens of kilometers long, across the eastern flank of the massif, encompassing the largest elevation gradient possible in the area, from the low-relief, high-elevation surface down to the bottom of the Rhône River valley, following the classical sampling strategy used for plateau margin studies (i.e., Persano et al., 2002). Analytical protocols for apatite AHe and AFT analysis adopted in this study followed Mahéo et al. (2013) and Balestrieri et al. (2011), respectively, and the corresponding procedures are presented in the Appendix.

Bedrock thermal histories were investigated by inverse and forward numerical modeling using QTQt and HeFTy software (Gallagher, 2012; Ketcham, 2005). QTQt software was used for inverse modeling, where observed AFT and AHe data such as spontaneous, induced, and horizontal confined tracks, apatite crystal dimension, amount of uranium and thorium, and measured age are inverted to find the temperature-time (T-t) paths compatible with the data. QTQt and HeFTy software was used in a forward approach, where different thermal histories are proposed following some specific erosional scenarios, and predicted data are directly compared to observed data.

The Ketcham et al. (2007) annealing model and the radiation damage accumulation and annealing model of Gautheron et al. (2009) were applied to the AFT and AHe data, respectively. In the case of a long stay in the He PRZ zone, special attention was paid to the impact of damage on He retention, as suggested by Ault et al. (2019).

Fourteen bedrock samples of granite and high-grade metamorphic rocks were collected along the eastern side of the Massif Central (Fig. 3). Single-crystal apatite AHe dating was performed on all samples, with three to five replicates resulting in 58 single-grain AHe ages (Table 1). Single-grain AHe variability within individual samples reflects variable grain size, damage content created during alpha decay, chemical composition, and complex diffusion processes (Shuster et al., 2006; Flowers et al., 2009; Gautheron et al., 2009, 2013; Recanati et al., 2017). AFT dating and track-length measurements were performed on five samples (Table 2; GSA Data Repository Fig. S1¹). In general, the AHe/AFT ages showed a positive correlation with elevation (Fig. 3), and the regression of the AHe age-elevation profile yielded an apparent exhumation rate of 0.04 km/m.y.

The northern profile was composed of four bedrock samples collected from elevations between 471 m and 1366 m asl, along the southern flank of St. Etienne valley. The valley is around 25 km wide, and the horizontal distance along our sampling profile was 12 km for a total elevation difference of 900 m (Fig. 3). Single-grain AHe ages span from 58 ± 5 to 131 ± 10 Ma. Five grains were dated for the lower three samples, while the uppermost sample yielded only one single-grain AHe age of 131 ± 11 Ma. Two AFT analyses from the highest and lowest samples showed ages of 93 ± 4 Ma (vo-14, 1136 m asl) with a long MTL of 13.9 ± 0.1 µm (100 measured tracks) and 75 ± 4 Ma (vo-16, 471 m asl) with an MTL of value of 12.7 ± 0.2 µm (51 measured tracks), respectively.

The southern profile was composed of seven bedrock samples from elevations of 195 m to 922 m, along a horizontal distance of ∼40 km, starting from the low-relief upland surface to the Rhône River, for a total difference in elevation of ∼700 m. Single-grain AHe ages span from 40 ± 3 to 133 ± 11 Ma. The three lowermost samples (vo-1, vo-2, vo-3) showed similar single-grain AHe ages with many young replicates (i.e., between 94 and 40 Ma). Two samples at 350 m (vo-4) and 220 m (vo-3) had AFT ages of 80 ± 4 Ma with MTL of 12.9 ± 0.2 µm (100 measured tracks) and 64 ± 4 Ma with MTL of 13.3 ± 0.2 µm (100 measured tracks), respectively. Sample vo-3, showing the youngest AFT age, failed the χ² test (P[χ²] < 5%) and is significantly dispersed (D = 19%; Table 2; Fig. S1). A single-grain age distribution failing the χ² test and showing such a high dispersion may be indicative of a complex thermal history (i.e., it resided for some time at partial annealing temperatures before ultimately cooling down during the final exhumation event that brought it to the surface). The young single-grain AFT ages corresponding to the second curve of the abanico plot in the radial plot (Fig. S1) indicate that, in the case of sample vo-3, the final exhumation event would have been in the Cenozoic.

Three bedrock samples were collected along a transect in between the two other profiles (vo-18, vo-19, vo-21), and they gave single-grain AHe ages spanning between 55 ± 4 and 129 ± 10 Ma. Sample vo-21 is the lowermost sample collected in this study, and it showed a high dispersion in single-grain AHe ages. The corresponding AFT age is 74 ± 4 Ma with MTL of 13.3 ± 0.1 µm (100 measured tracks).

The spread of single-grain AHe ages reflects different influences on the effective uranium content (eU, defined as U + 0.235 × Th), the grain size (sphere equivalent radius, Rs), and the thermal history (e.g., Gautheron et al., 2009). The apatite eU content was found to be highly variable between the different crystals (Table 1), but within individual samples, it was often positively correlated with the AHe ages, suggesting an important radiation damage control on AHe age dispersion (Fig. S1). This eU-AHe age trend is similar to the one obtained by Gautheron et al. (2009) for samples collected some tens of kilometers southward, where the AFT-AHe age overlap has been proposed to reflect both the level of radiation damage and the thermal history. In case of reheating or a long stay in the He PRZ, a high level of α-recoil damage can move the He PRZ to temperatures as high as 120 °C, resulting in a similar closure temperature for both the AHe and AFT systems (Gautheron et al., 2009).

To evaluate the thermal evolution of individual samples, taking into account both the influence of the α-recoil damage and grain size for the AHe system as well as AFT data, we performed one-dimensional thermal modeling using the QTQt software (Gallagher, 2012). Figure 4 shows the results of the modeling for the five samples from which AFT data and AHe ages were available. For inverse QTQt modeling, we used the radiation damage model for He diffusion (Gautheron et al., 2009) with the timing of the deposition of Triassic and Jurassic limestones as the only external stratigraphic constraint.

The thermal history of vo-14 (the uppermost sample) was modeled in two different ways, with and without the stratigraphic constraint (red lines in Fig. 4B). We note that the Mesozoic stratigraphic constraint does not modify the post–100 Ma thermal histories, suggesting that Mesozoic burial under Triassic and Jurassic sediments is not constrained by the thermochronological data for the uppermost sample.

For the low-elevation samples, the general thermal history showed a similar trend, with faster cooling phases from 120 to 100 Ma and from 70 to 50 Ma, followed by a less well-defined evolution suggesting a long permanence around 50–40 °C and then a final cooling phase until present day. It is noteworthy that the AHe data for the lowermost samples (vo-3 and vo-21; Figs. 4C and 4D) allow the thermal history to be substantially constrained, resulting in very similar outcomes between the two samples. In the case of no reliable AHe data, as for sample vo-4 (too much dispersion in the single-grain He ages), the Cenozoic thermal history cannot be precisely defined (Fig. 4E).

The geodynamic episodes and mechanisms that shaped the topography of the Massif Central are unclear because of the limited available time constraints for tectonic deformation events and geomorphological markers. The margin of the massif records a long morphogenesis history (see Séranne et al., 2002, and references therein), evolving under predominantly slow erosion rates. This is also supported by our new thermochronological data, which preserve the record of the Cretaceous exhumation but also reveal possible evidence for a second Cenozoic pulse. Cenozoic uplift and consequent erosion (and landscape rejuvenation) have been postulated by many authors (Barbarand et al., 2001), although quantification and timing have remained unconstrained. In the next section, we will trace the possible Cenozoic uplift signal in our thermochronological data and in published studies (Gautheron et al., 2009), and we will try to quantify the Cenozoic erosion from this data set. We remark that the slow erosion rates and limited amount of erosion prevent a straightforward interpretation of the thermochronological data, which are not univocal.

Comparing our new data set with the published data set of Barbarand et al. (2001) and Gautheron et al. (2009), we found some common features that we consider essential for the interpretation of the eastern margin evolution: (1) Uppermost samples show Cretaceous AFT and AHe ages and long MTLs (>13 µm); (2) a complex Cenozoic thermal history is suggested by the inverse modeling for the lowermost samples (in particular, for our samples vo-16 and vo-21; Figs. 4B and 4C); (3) MTL-elevation relationships show a complex and nonlinear trend (i.e., a general decrease of MTL from high elevation/old age toward intermediate elevations and then a slight increase for the lowermost and youngest samples; Fig. 5); and (4) large differences in AFT and AHe ages (∼40 m.y.) occur between the uppermost and lowermost samples despite the limited difference in elevation (∼1000 m).

The AFT age and MTL of our uppermost sample (vo-14; 93 ± 8 Ma and MTL of 13.9 ± 0.1 μm), located at 1136 m asl, are consistent with sample no. 20 (105 ± 6 Ma, MTL of 13.76 μm at 1500 m) from the Cévennes and with sample no. 8 (119 ± 4 Ma, MTL of 13.71 μm at 1185 m) from the Ardèche reported in Barbarand et al. (2001). The single-grain AHe ages of our sample vo-14 and sample no. 20 from Cévennes overlap too, being between 88 and 93 Ma, respectively. These new ages support the interpretation of Barbarand et al. (2001) and Peyaud et al. (2005) for a relatively fast and widely diffuse Cretaceous erosion phase, possibly induced by regional uplift associated with the Pyrenean rifting event. Initiation of this phase could have been around 110 ± 10 Ma, as derived by the thermal modeling of our data, which is consistent with Barbarand et al. (2001). The observed regional variations could reflect along-margin changes in tectonic evolution, which will be discussed below.

The lowermost samples from our study area show a complex Cenozoic thermal history, and the inverse modeling did not allow a complete fitting between observed and modeled data (Fig. 4). We used the forward modeling approach (QTQt with radiation damage diffusion model of Gautheron et al., 2009) to explore the three data sets (AFT, MTL, and AHe) and predictions in response to three possible thermal histories corresponding to different scenarios (margin erosion style and massif evolution): (1) constant cooling from Cretaceous to present, representative of a long-lasting topographic degradation (associated with erosion-related isostatic uplift; Figs. 5A and 5B); (2) a Cretaceous cooling phase ending in the mid-Cenozoic and representative of a single Cretaceous tectonic phase (Figs. 5C and 5D); and (3) a Cretaceous cooling phase followed by a phase of thermal stability around 40 °C and renewed (lower amplitude) cooling during the late Cenozoic, representative of two main tectonic phases, the first in the Cretaceous and a minor one in the late Cenozoic (Figs. 5E and 5F).

The three lowermost samples of our new data set have the highest probability of recording a possible Cenozoic phase, and thus they were modeled for the two AFT/AHe thermochronometers. Figures 5A and 5B show that for constant cooling since the Cretaceous, the predicted AFT age and MTL are similar to the observed ones, while the predicted AHe ages are too young for all three samples. The thermal history with a single Cretaceous event (Figs. 5C and 5D; similar to that proposed for the uppermost samples) tends to make the AHe ages older and similar to the observed ones, but it also predicts MTLs that are too long, making this scenario unlikely for the lowermost samples. The cooling history with both Cretaceous and Cenozoic cooling events allows us to reproduce the AFT age and MTL well, while it also tends to minimize the difference between the predicted and observed AHe ages, at least for the vo-3 and vo-21 samples (Figs. 5E and 5F). In general, the third scenario seems to best reproduce our overall data set.

The relationship of MTL with elevation and AFT ages is similar in all the three regions of Cévennes, Ardèche, and our study area (called the NE margin); the MTLs tend to decrease from the highest elevation to ∼600–400 m, and then they slightly increase from 400 m to the bottom of the profile (Fig. 6). This trend is highlighted by the second-order polynomial regression curves (calculated with age vs. elevation and then vertically rotated). The same trend is observed in the MTL versus AFT ages (Fig. 6D), where the Cévennes (and partially for the Ardèche and our data too) age patterns show a typical concave-up “boomerang” shape, which is commonly observed in passive continental margins (Green, 1986; Gallagher and Brown, 1999). Although the overall geodynamic evolution of the eastern Massif Central is very different from a continental margin, we note the similarity in terms of erosional style represented by a progressive erosive wave from a base level toward the interior that could produce comparable tectonic and geomorphologic evolution.

In the classical interpretation of these boomerang shapes for continental margins (Gallagher and Brown, 1999; Wildman et al., 2019), the upper samples preserve the erosional history of the premodern margin formation, in our case, the Cretaceous Pyrenean rifting (blue area in Figs. 6A–6C), while the lower samples record the evolution of the modern escarpment (green area in Figs. 6A–6C). In our case, the ages of the low-elevation samples do not precisely constrain the onset of erosion, because exhumation was not enough to expose samples from below the PAZ, as confirmed by MTLs that do not exceed 14 μm. From these ages, we can only propose a minimum age for the late-stage erosional history, which should be younger than 45 Ma, i.e., the youngest AFT age (45 ± 2 Ma) accompanied by an MTL of 12.91 µm (Barbarand et al., 2001).

In the age-elevation relationships, the AFT ages and the youngest AHe single-grain ages of each sample show a positive correlation with elevation (Fig. 3). In the case of Massif Central topography formed in the Cretaceous, the observed large difference in age for such a small difference in elevation is not expected because of the isotherm’s shape at depth: Isotherms are influenced by topographic wavelength, relief amplitude, and by isotherm depth (Stüwe et al., 1994; Braun, 2002; Foeken et al., 2007), so that below wide valleys (wavelengths >20–40 km), shallower isotherms tend to be parallel to the surface. If the topography of the Massif Central formed in Cretaceous to Eocene times, the exhumation of the investigated bedrock samples should have occurred under a topography similar to the present-day topography. Therefore, a sampling profile along a valley flank of more than 10–20 km long with less than 1 km difference in elevation is expected to yield a poor correlation between age and elevation (Braun, 2002).

The observed positive correlations in the AFT and AHe age-elevation relationships from the Massif Central suggest that relief was formed after these samples cooled below their closure temperature, providing further evidence in support of a Cenozoic exhumation event.

While in the previous section, we singularly tested the thermal history of the three lowermost samples, we also wanted to test if a Cenozoic exhumation phase and progressive erosion moving westward from the Rhone River to the interior would be consistent with the regional patterns of the thermochronological data. In this forward modeling approach (using HeFTy; Ketcham, 2005), we iteratively defined T-t histories to reproduce the observed AFT and AHe ages and MTL distributions along the NE margin area E-W profile. Then, we did the same for the published data along the Cévennes profile, but, there, we did not include the AHe ages because multigrain samples were used in the study by Gautheron et al. (2009), and so the comparison with single-grain predicted age is too uncertain.

In Figure 7A, we present a swath profile running perpendicular to the NE margin. We integrated the new AFT and AHe data with two published samples (no. 1 and no. 2 of Barbarand et al., 2001). We used the forward modeling to find a T-t path producing data as similar as possible to the obtained AFT data and AHe ages. We report the forward modeling path with the inversion modeling obtained using QTQt software (Fig. 7B) for comparison.

In the NE margin, the lowest samples vo-21 and vo-3 showed MTL longer than the lowermost samples of the Cévennes and Ardèche regions. The T-t path constrained for samples vo-21 and vo-3 is a burial (∼110 °C) during the late Mesozoic followed by cooling to shallow crustal levels corresponding to ∼40 °C, where the sample resided for ∼20 m.y. before being exhumed since 40–30 Ma. Interestingly, this proposed exhumation history is consistent with the thermal histories from QTQt modeling, which also relied also on single-grain AHe ages (colored paths in Fig. 7). Sample vo-4 yielded a shorter MTL, consistent with a longer stay at temperatures ∼50 °C and a younger cooling event since ca. 30 Ma. Although the two samples (vo-3 and vo-4) are quite close, the distance of ∼5 km is considered enough to explain the observed difference in MTL as being the result of progressive erosion propagation that had reached the most internal sample a few million years later. Considering the distance of ∼5 km between the two samples, we can suggest a first-order propagation rate of the erosion wave between ∼0.1 and 1 km/m.y., which is consistent with the escarpment retreat rate found in other continental margins (Braun, 2018; Wildman et al., 2019).

Sample vo-16 showed data similar to sample vo-4 and occupied almost the same location, even if it was located further north, where the margin is steeper, and its projection is reported on our profile. Thus, the envisaged T-t path is composed of Cretaceous cooling to temperatures of 30–40 °C and, after some million years of stagnation, a renewed phase of cooling starting at ca. 20 Ma. Sample vo-14 is aligned with sample vo-16, but it is almost on top of the massif at higher altitude, and its data are well reproduced by a Cretaceous cooling phase bringing it almost to surface temperatures.

Most internal samples (no. 2 and no. 1 of Barbarand et al., 2001) were not completely reset during the Cretaceous burial and preserve the inheritance of the old exhumation histories. Moreover, because sample no. 2 is a Triassic sandstone, it further supports limited postsedimentation erosion. In general, their short MTLs are not consistent with a Cenozoic event suggesting slow cooling. A very slight change in cooling rate at ca. 30 Ma can be observed in the modeling of sample no. 1 only (Fig. 7).

Next, we applied the forward modeling to the Cévennes profile (Barbarand et al., 2001) for comparison (Fig. 8). We chose seven samples from Barbarand et al. (2001) (no. 20, 21, 22, 25, 26, 27, and 28) showing the AFT and track-length distributions considered as representative of the thermal histories previously discussed. Thermal histories of low-elevation samples no. 27 and no. 21 (Fig. 8), which gave AFT age of 51 ± 2 and 45 ± 2 Ma and MTLs of 12.30 and 12.91 µm, respectively, are consistent with a limited amount of burial during the Cretaceous (up to ∼96 °C), stable or slow exhumation since ca. 90 Ma, and an increase in cooling rate starting at ca. 40 Ma from a temperature of ∼80 °C. For these samples, the Cretaceous cooling phase left only a weak signal, likely because, at that time, they were buried at great depth, and they did not record this erosion event. Samples no. 22 and no. 26 of Barbarand et al. (2001) were collected from further into the interior and should have recorded the passage of the Cenozoic erosional wave some million years later due to the westward propagating erosion from the base level to the margin top. Samples no. 22 and no. 26 seem very close in Figure 8 because they are projected on the swath profile, but they are ∼10 km from each other. Shorter MTLs for samples no. 22 and no. 26 are consistent with the model of a westward propagating erosion wave, a longer residence time in the PAZ, and thus a more recent onset of exhumation. The onset of exhumation and Cenozoic increase in cooling rate started when samples no. 22 and no. 26 were at ∼70 °C, and, consequently, they yielded a mixed age with some inheritance from the pre-Cenozoic history. Samples no. 25 and no. 28 recorded old cooling periods (120–80 Ma) due to the effects of the Cretaceous erosion event, even though it was not of sufficient magnitude to bring the samples to the surface, such as suggested by short MTL (very short for no. 28). Their last cooling phase (from ∼60 °C to surface temperatures) was induced by the Cenozoic event. For these samples, the onset of this last cooling phase was even younger due to the progression of the erosion wave.

The uppermost sample, no. 20 (Barbarand et al., 2001), showed a fast cooling rate starting ca. 120 Ma, while the pre-Cretaceous history was totally reset and remains unconstrained (dashed lines in the modeling of sample no. 20). Sample no. 20 does not record any evidence for Cenozoic erosion because, after the Cretaceous phase, it was already near surface temperatures.

The Cévennes, Ardèche, and NE margin (our study) areas were part of the western margin of the Southeast France Basin during the Cretaceous. This sedimentary basin developed on European crust since the Triassic in response to Alpine Tethys rifting. Barbarand et al. (2001) interpreted the mid-Cretaceous AFT cooling age pattern as a general erosional event produced by an uplift commonly referred to as the “Durance uplift” event, which involved the Southeast France Basin and probably was associated with the opening of the Bay of Biscay. Here, we propose that mid-Cretaceous thermochronological ages not only recorded the “Durance uplift” event, but they might also reflect the specific evolution of the Southeast France Basin.

The Southeast France Basin was initially NE-SW elongated in the Triassic to Middle Jurassic, being parallel to the northern Tethyan paleomargin; in the Late Jurassic, the structural trend became E-W. The late Early Cretaceous was a phase of deep reorganization in the Southeast France Basin geometry, associated with general N-S–directed extensional tectonics that led to the formation of the coeval “isthmus of Durancien” and two separate basins, the Vocontian basin in the north and the South Provence basin in the south. It is possible that the erosional event recorded in the thermochronological data was induced by the margin uplift at the transition between a subsiding and extending continental domain in the east (the Vocontian Basin) and the Massif Central domain to the west.

During the Barremian (129–125 Ma), the westward portion of the Vocontian trough was close to the Ardèche area, while during the Albian (ca. 110–100 Ma), the Vocontian trough expanded southward, toward the Cévennes, and northward toward the NE margin (our study area; Curnelle and Dubois, 1986). The younger AFT/AHe ages found in the uppermost samples in the Cévennes and in our study area (ca. 120–90 Ma) with respect to the older Ardèche samples (ca. 130–110 Ma) might reflect the nonsynchronous uplift of the margin that progressively followed the basin subsidence.

In such a scenario, the NE-SW–oriented margin of the Massif Central was oblique with respect to the general N-S–directed sense of extension, producing a possible transtensional tectonic regime. The thickness of removed crust (several kilometers) along the margin is indeed compatible with a mechanism of erosion related to a transtensional uplifted rift shoulder similar to that observed in modern continental margins. Such a geodynamic context is similar to propositions for oblique continental margins such as the Transantarctic Mountains (Wilson, 1995).

The importance of the Cretaceous extensional phase, associated with the Bay of Biscay–Pyrenean rifting, is supported by a large body of geological evidence in the Southeast France Basin and also recently found in the western Alps (Tavani et al., 2018, and references therein).

The hypothesis of a Cenozoic age for the eastern Massif Central topography is supported by the thermal and topographic evolution derived from our data set and published data. These data sets and associated modeling outcomes are consistent with a regressive erosion wave moving from the topographic margin to the massif interior following a Cenozoic uplift phase.

The massif margin may be considered as a smaller-scale analogue to an elevated continental passive margin, with an erosional pattern due to a regressive erosion wave propagating from the base level to the top following a Cenozoic uplift phase (see Wildman et al., 2019). The base level is represented by the Rhône River, which contributed to maintain the main valley floor at the same elevation during the Cenozoic. The limited amount of erosion, coupled with the occurrence of Cretaceous deposits on top of the massif and in the Rhône river valley floor and no evidence for faulting, suggests long-wavelength flexure of the lithosphere, which produced a margin topography characterized by a broad monocline with a very low gradient. In the case of the Massif Central, no Cenozoic extensional rifting domain is invoked to produce the observed uplift, but such an erosional pattern is consistent with long-wavelength topographic growth induced by mantle upwelling, providing a remarkable consistency with the deep and localized present-day high-temperature anomaly imaged by seismic tomography (Granet et al., 1995; Faccenna et al., 2014) and the long-term surface evolution. Volcanism is further evidence for mantle involvement during topographic growth, with the first activity recorded ca. 30 Ma and an effusive acme since 15 Ma (Michon and Merle, 2001).

The onset time for uplift is not well constrained by our thermochronology data. The Cenozoic event recognized in the MTL regional distribution and supported by forward modeling (Figs. 5 and 6) should have occurred after ca. 45 Ma (i.e., the AFT age of the youngest sample in the Cévennes area; Barbarand et al., 2001), although it remains difficult to precisely constrain because of limited erosion and dispersed single-grain AHe ages. Lutetian (48–40 Ma) marine deposits have been described in several basins located in the middle of the massif, such as the Puy Basin (Fig. 2; Séranne et al., 2002; Michon, 2000, and references therein), suggesting that margin uplift did not occur before ca. 40 Ma (Fig. 9A).

Between the late Eocene and Oligocene, some portions of the Massif Central were involved in two regional extensional events: European rifting and the opening of the Gulf of Lion. The European rifting was associated with localized crustal extension within individual basins, such as the Bresse basin, while a coeval uplift and erosional event at regional scale has been never described. The amounts of uplift and volcanic activity in the Massif Central are much higher than in other analogue massifs (e.g., Vogelsberg mountains, Rhenish and Bohemian Massifs), suggesting the contribution of other processes. The opening of the Gulf of Lion (30–16 Ma; Fig. 9C) corresponds to an important tectonic phase that reactivated Paleozoic structures such as the Nimes fault and in general deeply reshaped the western Mediterranean plate organization (Séranne et al., 1995). The onset of margin uplift could have started at the beginning of this reorganization event.

The Massif Central plume in the early Miocene probably interacted with the Ionian subduction system, which contributed to continuously nourish the mantle upwelling and flow. The contribution of Ionian subduction to Massif Central volcanism has been already proposed (Barruol and Granet, 2002; Faccenna et al., 2010) and also supported by seismic anisotropy (Faccenna et al., 2014; Salimbeni et al., 2018), which imaged the NW-SE–directed mantle counterflow induced by Ionian slab retreat during the opening of the Gulf of Lion (Fig. 8C). The interaction of the Massif Central mantle plume with the slab retreat–induced counterflow could have caused the greater volcanism, the higher total elevation, and the long-term slow surface uplift evolution of the Massif Central.

The amounts of exhumation and margin uplift remain uncertain. If we consider the thermal modeling consistent with progressive westward erosion of the margin (Figs. 7 and 8), we can observe that samples cooled by ∼40 °C during the last 30–40 m.y., corresponding roughly to the removal of ∼1 km of crustal thickness (assuming a mean annual surface temperature of 10 °C and a geothermal gradient of 30°/km). This results in an exhumation rate ranging from 0.02 to 0.03 mm/yr (Fig. 9C), which is comparable with modern basin-averaged cosmogenic denudation rates reported in the area (∼0.04 mm/yr; Olivetti et al., 2016), and which suggests that a large portion of the slowly eroding landscape underwent a constant process of erosion during the last 40 m.y. Other portions of the landscape (mainly in the lowermost portion of the catchments) are eroding slightly faster (cosmogenic denudation rate of 0.07 mm/yr), suggesting that uplift rate might have further increased, and the uplift is probably still active today. This possible recent or ongoing uplift is responsible for the elevation of Pliocene marine-continental transition deposits found along the Rhône valley, at 200 m asl (Peage de Roussillon village in Fig. 2; Aguilar et al., 1989).

Our new thermochronological data from the eastern margin of the French Massif Central, combined with published data, enabled us to reconstruct the long-term evolution of long-wavelength topography. A two-step exhumation history has been recorded in the regional data trends: (1) a first mid-Cretaceous exhumation phase, and (2) a second, newly recognized, Cenozoic erosional phase affecting the newly formed margin with onset time likely after the Eocene. For this Cenozoic phase, we envisage a domal uplift of the area at least partly due to mantle upwelling, with the formation of a long-wavelength flexure at the margin and subsequent erosion of the scarp. The total amount of erosion is limited, such that AHe and AFT ages exposed by the erosional wave do not quantitatively constrain the onset of the phase but rather yield mixed ages. On the contrary, forward and inverse thermal modeling of AFT and AHe ages plus track-length distribution permitted us to unravel the signal of a post–50 Ma event of erosion induced by a progressive westward regressive erosion starting at the eastern margin of the massif.

The research of and postdoctoral allocation to Olivetti were funded by Electricité de France (EDF) through the SIGMA research program (SeIsmic Ground Motion Assessment). This work is also a contribution of the Ecosystèmes Continentaux et Risques Environnementaux (ECCOREV) research federation and Labex OT-Med (ANR-11-LABX-0061) funded by the French Government “Investissements d’Avenir” program of the French National Research Agency (ANR) through the A*MIDEX project (ANR-11-IDEX-0001–02). Valla acknowledges support from Swiss National Science Foundation (SNSF) grant PP00P2_170559. We thank Science Editor Damian Nance and three anonymous reviewers for constructive criticism and suggestions.

The apatite grains were recovered from the collected bedrock samples following standard crushing, sieving, washing, magnetic, and heavy liquid separation.

The AHe analyses were carried out at the Paris-Sud University. Euhedral apatites were picked using a cross-polarized binocular microscope. Most grains had a minimum diameter of 90 μm and were inclusion free to avoid effects of He-implanting from inclusions or excess loss of He during decay due to a large surface/volume ratio (Farley, 2000). The grain dimensions were measured for calculation of the alpha-ejection (Ft) correction factor after Farley et al. (1996), and single grains were packed in Nb tubes for U-Th/He measurement. For each sample, up to five aliquots were prepared for analysis in order to ensure sample age reproducibility. The concentration of ⁴He was determined by the ³He isotope dilution and measurement of the ⁴He/³He ratio with a quadrupole mass spectrometer. Apatite samples were heated for 5 min at 11 Amps with a 960 nm diode laser for degassing. Each sample was reheated and measured to ensure that all gas was extracted in the first run. U and Th concentrations were obtained by isotope dilution using an inductively coupled plasma–mass spectrometer.

Apatite grains were mounted in epoxy resin, ground, and polished to expose internal mineral surfaces. Etching with 5 N HNO₃ at room temperature for 20 s revealed spontaneous fission tracks intersecting the apatite surface. Samples were covered with a uranium-free muscovite external detector and irradiated with thermal neutrons at the LENA (Laboratorio Energia Nucleare Applicata) Triga Mark II reactor of Pavia University, Italy. Induced fission tracks in the external detector were revealed by etching the mounts in 40% HF at room temperature for 40 min. The fission tracks were counted by the first author under a nominal magnification of 1250× on a Zeiss Axioskop equipped with a Kinetek automatic stage at l’Istituto di Geoscienze e Georisorse (IGG) fission-track laboratory, an institute of the Consiglio Nazionale delle Ricerche (CNR). The Trackkey 4.2 Program was used for all AFT age calculations procedures (Dunkl, 2002). A chi-square (χ²) test was carried out on the AFT single-grain ages in order to test the homogeneity of data (Galbraith, 1981). The probability of (χ²) was calculated for each sample; if P(χ²) > 5%, then the sample was assumed to be homogeneous (Galbraith and Laslett, 1993).