Abstract

Resolving the timing of initiation and propagation of continental accretion associated with increasing topography and exhumation is a genuinely challenging task using low-temperature thermochronology. We present an integrated thermo-mechanical and low-temperature thermochronology modelling study of tectonically-inverted hyperextended rift systems. Model low-temperature thermochronology data sets for apatite (U-Th)/He, apatite fission-track, zircon (U-Th)/He and zircon fission-track systems, which are four widely used thermochronometric systems in orogenic settings, are generated from fourteen locations across a model collisional, doubly-vergent orogen. Our approach allows prediction of specific, distinct low-temperature thermochronology signatures for each domain (proximal, necking, hyperextended, exhumed mantle) of the two rifted margins that, in turn, enable deciphering which parts of the margins are involved in orogenic wedge development. Our results show that a combination of zircon (U-Th)/He and apatite fission-track data allows diagnostic investigation of model orogen tectonics and offers the most valuable source of thermochronological information for the reconstruction of the crustal architecture of the model inverted rifted margins. The two thermochronometric systems have actually very close and wide closure windows, allowing to study orogenic processes over a larger temperature range, and therefore over a longer period of time. Comparison of model data for inverted rifted margins with model data for non-inverted, purely thermally-relaxed rifted margins enables assessing the actual contribution of tectonic inversion with respect to thermal relaxation. We apply this approach to one of the best-documented natural examples of inverted rift systems, the Pyrenees. Similarities between our thermochronometric modelling results and published low-temperature thermochronology data from the Pyrenees provide new insights into the evolution of the range from rifting to collision. In particular, they suggest that the core of the Pyrenean orogen, the Axial Zone, consists of the inverted lower plate necking and hyperextended domains while the Pyrenean retrowedge fold-and-thrust belt, the North Pyrenean Zone, represents the inverted upper plate distal rifted margin (exhumed mantle, hyperextended and necking domains). This is in good agreement with previous, independent reconstructions from literature, showing the power that our integrated study offers in identifying processes involved in orogenesis, especially early inversion, as well as in predicting which domains of rifted margins are accreted during mountain building.

Résumé

Identifier au moyen de la thermochronologie basse température le moment à partir duquel l’accrétion continentale s’amorce et se propage est une tâche difficile. Dans cette étude, nous couplons modélisation thermomécanique de systèmes hyper-amincis inversés et thermochronologie basse température pour prédire la signature thermochronologique de quatre systèmes basse température communément utilisés en domaine orogénique ((U-Th)/He sur apatite, traces de fission sur apatite, (U-Th)/He sur zircon et traces de fission sur zircon). Notre approche de modélisation permet de prédire les signatures thermochronologiques basse température propres à chaque domaine de marge du système extensif (domaine proximal, zone de neck, domaine hyper-aminci et domaine de manteau exhumé). Ces signatures permettent en retour d’identifier les domaines de marge impliqués dans l’accrétion continentale qui suit l’extension. Nos résultats montrent que la combinaison des signatures thermochronologiques des systèmes (U-Th)/He sur zircon et traces de fission sur apatite renseigne sur les processus mis en jeu lors de la formation d’un orogène collisionnel à double vergence, ainsi que sur l’architecture des marges inversées durant l’orogénèse. Les deux systèmes thermochronométriques possèdent en effet des gammes de température de fermeture voisines et relativement larges, permettant d’étudier les processus orogéniques sur un intervalle de températures plus grand, et donc sur une période de temps plus longue. La comparaison des signatures thermochronologiques prédites dans un modèle thermomécanique d’inversion avec celles prédites dans un modèle de relaxation thermique, sans convergence postérieure à l’extension, permet d’apprécier la nécessaire contribution des processus tectoniques d’inversion par rapport au processus de relaxation thermique postrift. Nous appliquons notre approche à l’un des exemples naturels les mieux documentés de systèmes de rift inversé, les Pyrénées. Les similitudes qui existent entre nos prédictions thermochronométriques et les données de thermochronologie basse température publiées pour les Pyrénées fournissent de nouvelles informations sur la formation et l’évolution de l’orogène pyrénéen. En particulier, ces similitudes suggèrent que la zone interne de déformation pyrénéenne, la Zone Axial, est constituée des domaines de neck et hyper-aminci de la plaque inférieure, inversés durant l’accrétion, tandis que le prisme orogénique qui se développe sur la plaque européenne, la Zone Nord-pyrénéenne, représente la marge distale inversée de la plaque supérieure. Ce résultat est en accord avec les précédentes reconstructions publiées de la dynamique pyrénéenne, démontrant ainsi le potentiel de notre approche dans l’identification des processus impliqués dans l’orogenèse, notamment l’inversion précoce, ainsi que dans la prédiction du type de domaine de marge accrété lors de la formation des chaînes de montagnes.

Introduction

Thermally-coupled dynamical modelling of low-temperature thermochronology data in orogenic wedges has proved to be useful to quantify slip rates, dip and vergence of thrusting, to define tectonic scenarios and processes involved in continental accretion, and to examine the role of external climate forcing into mountain building in a variety of settings (Taiwan, the Himalayas, New Zealand, the Pyrenees, the Andes; e.g.Batt et al., 2001; Willett et al., 2003; Herman et al., 2010; Braun et al., 2012; Erdös et al., 2014; McQuarrie and Ehlers, 2015). This approach has been preferentially focused on rapidly exhuming orogenic systems (e.g. Taiwan, the Himalayas, New Zealand) where efficient, coupled tectonic-geomorphic processes provide measurable changes in low-temperature thermochronology signals. Yet, these processes are expected to be less significant during incipient collision, when tectonic inversion of thin distal rifted margins occurs under water with little to no erosion. In this case, cooling recorded by low-temperature thermochronology would rather reflect deep tectono-thermal processes. It is suggested, for instance, that the downward deflection of isotherms related to fast underthrusting may alone explain the early, syn-convergence cooling event documented in the Taiwan accretionary prism (Mesalles et al., 2014). Alternatively, in low-convergence orogenic systems like the Pyrenees, both thrusting and underplating processes are shown to be able to generate rapid cooling during inversion of the distal margins of narrow rifts (Ternois et al., 2019a) while sedimentary burial over large hyperextended domains can maintain temperatures above the sensitivity limit of low-temperature thermochronometers and delay cooling record until later collision of the two proximal margins (Vacherat et al., 2014, 2016). Eventually, where the orogenic prism preserves the thermal history of the rifted margins, low-temperature thermochronology data may reveal cooling due to postrift relaxation rather than extensional unroofing (Malusà et al., 2016; Rat et al., 2019). Taken together, these independent propositions therefore show the non-unique significance of low-temperature thermochronology signals in orogens. This is observed not only from one collisional setting to another (e.g. Taiwan compared to the Pyrenees) but also from one domain of a rifted margin to another (e.g. the distal domain compared to the hyperextended domain in the Pyrenees). Firstly, this questions the potential of low-temperature thermochronology to faithfully record cooling related to the onset of orogenesis, depending on the crustal architecture of the accreted rifted margins. Secondly, this raises the need for thermal models of mountain building that, unlike existing kinematic models, explicitly and self-consistently integrate the mechanical evolution of collisional orogens from rifting to collision. Not only this would allow full exploitation of the constraints provided by low-temperature thermochronology, but this would also enable evaluation of the respective contribution of crustal and surface processes during the early phases of continental accretion.

Here, we take advantage of the recent numerical modelling results of rift inversion, designed to examine the impact of variable rift maturity, especially the effect of residual rift thermal anomaly, on the structure and strain distribution in orogenic belts by Jourdon et al. (2019). In order to constrain the tectonic and thermal evolution of domains of rifted margins (proximal, necking, hyperextended, exhumed mantle) during rifting, early convergence and subsequent collision, we extract from the numerical experiment thermal histories that are representative of the spatial and temporal evolution of upper crustal basement particles in each domain (Model A; Figs. 1A and SI01). To gain insights into the role of postrift thermal relaxation we also run a new numerical experiment with the same characteristics as in Jourdon et al. (2019) but without inversion (Model B; Fig. 2). We then use individual particles’ thermal histories returned by the two thermo-mechanical models as input for forward modelling of low-temperature thermochronology data. Our integrated thermo-mechanical and low-temperature thermochronology modelling study allows evaluating the ability of different, widely used thermochronometric systems to decipher key orogenic processes and to unravel the role of rift inheritance into early orogenic processes. Our approach is ultimately used through a case study of the Pyrenees to help predicting which parts of the rifted margins (proximal, necking, hyperextended, exhumed mantle domain) are inverted and exhumed across the orogen based on low-temperature thermochronology signals.

Methods

We describe below our approach for modelling the coupled physical and rheological evolution of tectonically-inverted hyperextended rifts. We select material points at the top of the model basement at the end of the experiment. Depth and temperature of the selected points are extracted backward in time for each model time step. As modelling output, thermal histories of the material points, referred thereafter as “time-temperature paths”, thus take into account evolution of the geothermal gradient through time and space.

Thermo-mechanical numerical modelling

The thermo-mechanical models presented in this study (Model A in Figs. 1A and SI01; and Model B in Fig. 2) are performed with pTatin2D (May et al., 2014, 2015). The code solves the conservation of momentum for an incompressible fluid as follows:
∇⋅(2ηε˙)P=ρg,
graphic
(1)
where η is the viscosity of the fluid, P is the thermodynamic pressure, ρ is the density of the fluid, g is the vector of the gravity acceleration and ε˙forumla is the strain rate tensor computed as:
ε˙=12(v+vT),
graphic
(2)
where v is the velocity vector of the fluid.
The mass conservation is computed for an incompressible fluid taking into account the Boussinesq approximation:
.v=0.
graphic
(3)
The temperature of the domain is computed with the time-dependent advection-diffusion energy conservation equation such as:
Tt+v.T=.(κT)+HρCp,
graphic
(4)
where T is the temperature of the fluid, t is the time, v is the velocity vector of the fluid, κ is the thermal conductivity, ρ is the density of the fluid, Cp is the thermal capacity and H is the sum of the different heat sources.
In the simulations, the computed heat sources are the radiogenic heat production computed accordingly to Turcotte and Schubert (2002):
Hprod=H0exp(yyp),
graphic
(5)
where H0 is the heat production at the surface, y is the depth and yp is the depth of the radiogenic layer.
In addition, the heat production due to mechanical energy dissipation is computed as:
Hshear=2ηε˙IIε˙II,
graphic
(6)
where ε˙IIforumla is the second invariant of the strain rate tensor computed as:
ε˙II=12ε˙ijε˙ij.
graphic
(7)
To simulate the mechanical behavior of the lithosphere at geological timescales, the brittle parts of the lithosphere deform accordingly to the pseudo-plastic Drucker–Prager yield criterion:
ηp=Ccos(φ)+Psin(φ)ε˙II,
graphic
(8)
where C is the cohesion of the material, P is the pressure and φforumla is the friction angle.
To consider small-scale processes contributing to the mechanical softening along faults plane, a linear friction decrease with increasing plastic strain is applied such as:
φ=φ0εpεminεmaxεmin(φ0φ),
graphic
(9)
where φforumla0 is the initial friction angle, φforumla is the final friction angle, εp is the plastic strain, εmin is the minimum amount of plastic strain before friction drop starts and εmax is the maximum amount of plastic strain above which no more friction decrease applies.
In addition, the ductile parts of the lithosphere are simulated with the non-linear viscous Arrhenius flow law for dislocation creep:
ηv=A1n(ε˙II)1n1exp(Q+PVnRT),
graphic
(10)
where A, n and Q are material-dependent variables, R is the gas constant and V is the activation volume.

In the models, the sediments and the upper and middle crust are simulated with a quartz rheology (Ranalli and Murphy, 1987), the lower crust is simulated with an anorthite rheology (Rybacki and Dresen, 2000) and the mantle is simulated with a dry olivine rheology (Hirth and Kohlstedt, 2003).

We consider first the Jourdon et al. (2019) model called “Run 2” as it bears obvious similarities of duration and rates of extension/compression with a natural inverted rift system, the Pyrenees. This run lasts 92.75 Myrs and involves a constant extension velocity of 5 mm/yr over 30 Myrs, a 2.5-Myr period of linearly decreasing extension velocity (5 to 0 mm/yr) followed by a 2.5-Myr period of linearly increasing convergence velocity (0 to 2.5 mm/yr) to simulate smooth plate boundary inversion, and 57.75 Myrs of convergence at a constant velocity of 2.5 mm/yr (Model A; Figs. 1A and SI01). The rheological layering assumed in the model setup permits crustal decoupling during rifting and collision in the middle-lower crust. Jourdon et al. (2019) showed that variation in the width of the exhumed mantle domain and in the thermal regime under which the convergence stage initiates leads to very different mountain belt evolution and final structure. The orogenic wedge produced in the “Run 2” model displays a retrowedge that is devoid of major tectonic overprints due to collision. The margin structure is therefore well preserved, albeit buried below thick syn-orogenic deposits (Fig. 1A). The mantle exhumed between the two plates during extension is partially accreted while its main part plays the role of a buttress permitting the accretion of the downgoing plate. The final orogenic belt is composed of stacked basement units originating from the necking and proximal domains that are detached in the middle-lower crust.

To study the effect of postrift thermal relaxation onto low-temperature thermochronology signatures, we also run a new, convergence-free model (Model B; Fig. 2) computed with the model setup of Jourdon et al. (2019) for the first 32.5 Myrs of rift evolution. The last 60.25 Myrs of the simulation are modified from the original setup, with thermal reequilibration replacing smooth plate boundary inversion and convergence. This model serves as a reference model of rift evolution in our study.

Low-temperature thermochronology data modelling

In Model A we select 14 time-temperature paths across the model orogen (Tab. SI01), each of them characterizing the time-temperature history of a numerical particle chosen not only to be lying above mid-crustal levels before rifting (depths < 10–15 km) but also to be reaching near surface conditions during inversion (Figs. 1 and SI01). Particles 1–3 and 13–14 are originally found in the slightly stretched proximal margins. Particles 4–5 and 11–12 are located in the necking domains, characterized by crustal thinning from about 30 to 10 km. Particles 6–8 and 10 are taken in the hyperextended domains, defined by crustal thicknesses below 10 km. Regions of exhumed mantle are ultimately represented by particle 9. The same particles are then chosen for Model B and corresponding time-temperature paths are extracted.

We transform the forward time in the thermo-mechanical models into backward time, so that the time displayed on time-temperature path diagrams corresponds to an age with respect to the end of the simulation computed as:
tbackward=tforwardtfinal,
graphic
(11)
where tbackward is the backward time at each time step, tforward is the forward time at each time step and tfinal is the forward time at the end of the simulation (92.75 Myr). This means that 92.75 Myr forward time becomes 0 Ma backward time, and 0 Myr forward time becomes 92.75 Ma backward time.

We then use paths from the two thermo-mechanical models as thermal history input for forward modelling of low-temperature thermochronology data using the software program QTQt (Gallagher, 2012; v5.7.1). QTQt comprises a “forward model”, which allows prediction of the expected data distribution for any given set or parameter values, i.e. any given thermal history. On a practical level, the “forward model” describes the evolution of a given instance of a thermochronometric system as a function of time and temperature assuming a particular starting arrangement and subsequent time-temperature history. These instances include concentration of daughter products that result from radioactive decay of parent isotopes (4He nuclei; fission-tracks) and alteration or loss of these daughter products (alpha ejection; fission-track annealing/length reduction). The forward modelling approach has become possible after time-dependent, temperature-sensitive processes were measured in the laboratory and their behaviour confidently extrapolated to geological time scales. The reader can refer to Ketcham (2005) for further details. To resolve the thermal history of upper crustal rocks in a temperature window that is typical of pre-magmatic orogenic processes (0–300 °C; e.g.Carrapa, 2010; Peyton and Carrapa, 2013, and references therein), we consider the single dating techniques that have proved to be powerful when reconstructing the history of most orogens: apatite (U-Th)/He (AHe; closure temperatures of ∼ 40–100 °C assuming a 10 °C/Myr cooling rate; Farley, 2000; Shuster et al., 2006; Flowers et al., 2009; Gautheron et al., 2009; Shuster and Farley, 2009), zircon (U-Th)/He (ZHe; closure temperatures of ∼ 140 to 200–250 °C; Reiners et al., 2002, 2004; Stockli, 2005; Tagami et al., 2003; Wolfe and Stockli, 2010; Guenthner et al., 2013), apatite fission-track (AFT; closure temperatures of ∼ 60–110 °C; e.g.Ketcham et al., 2007b), and zircon fission-track (ZFT; closure temperatures of ∼ 180 to 250–320 °C; e.g.Bernet, 2009) systems.

As it is reported in most geological settings, low-temperature single-system data sets come with intrasample age variability (e.g.Gautheron et al., 2009; Guenthner et al., 2014, 2015). It has been demonstrated that this results from the inability of a single set of kinematic parameters for He diffusion or fission-track annealing to describe the natural thermochronometric behavior of apatite (e.g.Shuster et al., 2006; Flowers et al., 2007; Ketcham et al., 2007a; Gautheron et al., 2009) and zircon (e.g.Reiners et al., 2002, 2004; Cherniak et al., 2009; Guenthner et al., 2013; Ketcham et al., 2013). To reproduce such variability in model single-system AHe, AFT, and ZHe data sets, we consider system-specific characteristics that are recognized as sources of intrasample grain age variation by the geo-thermochronology community (e.g.Barbarand et al., 2003; Flowers et al., 2009; Guenthner et al., 2013). These characteristics include intraparticle bulk-grain effective uranium concentration and grain size for the ZHe and AHe systems and intraparticle grain chlorine content for the AFT system. Modelling details are presented in the Supplementary Material (Tab. SI02). Regarding the age variability one might expect in ZFT data sets, natural sources of predicted age variation for the ZFT system are, as yet, little known. We therefore consider only parameterization variation in radiation damage annealing models for the ZFT system (see models by Tagami et al. (1998) and Yamada et al. (2007); Tab. SI02). As our considerations encompass ranges of values that are larger than usually measured for the different, potential sources of age variation (especially a chlorine content of 0–6 wt% for the AFT system), we argue that the predictions presented in Figures 3 and 4 in the next section describe the full range of data that could be obtained in a low-temperature thermochronology study.

Results

Figures 3A and 3B show model thermochronological ages (ZFT, AFT, ZHe, AHe) and fission-track length data generated for each time-temperature path from Model A (Figs. 1 and 3C). Figures 4A and 4B show thermochronometric data predicted for Model B with time-temperature paths plotted in Figure 4C. These results are reported in Table SI03 in the Supplementary Material. Age predictions are represented as a function of particle number (the original position of the particles across the rifted margins; Figs. 1A and 2) and predicted mean track lengths (MTL) are plotted against model AFT age using the particle colour scheme used in Figure 1.

Thermochronological age distributions in inverted rift systems (Model A)

Model ZFT age distribution: record of prerift and synrift thermal histories

Rifted margins are hereafter referred to as lower and upper plates on the basis of their position within the orogenic system at the end of convergence (Fig. 1). All domains of the inverted rift system show model ZFT ages (73 to ca. 93 Ma) older than the age of convergence onset in Model A (ca. 60 Ma; Fig. 3A and Tab. SI03). This indicates that these model ages are not fully reset by burial and any source of heating during convergence, and rocks largely preserve the thermochronologic fingerprint > 250 °C acquired before inversion. The older ZFT age predictions (ca. 93 Ma; numerical modelling starting time) are synchronous with the time of development of extensional brittle deformation in Model A and almost exclusively expected in inverted upper and lower plate proximal and necking domains, along with long zircon MTL. These model ages mark the beginning of modelling (numerical onset of extension) and can be considered as pre- to early synrift cooling ages. In contrast, as shown by the younger mean model ZFT ages they yield, particles 5 and 6 in the lower plate necking domain (87 and 91 Ma, respectively), particle 8 in the lower plate hyperextended domain (73 Ma), and particle 9 in the upper plate exhumed mantle domain (81 Ma), record partial reset during rifting, and possibly during postrift and later phases too. This also manifests through shorter predicted zircon MTL for particles 5, 6 and 9 (Fig. 3B and Tab. SI03), although we acknowledge that in practice these data are difficult, if not impossible, to obtain with low uncertainty in zircons and are consequently little used in recent low-temperature thermochronology studies.

During extension in Model A, particles 5 and 6, originally lying at shallower crustal levels (2–3 km), record moderate heating followed by cooling down to temperatures of 150–50 °C (Fig. 3C) due to particle exhumation in the footwall of a normal fault (Fig. 1). Around 10 Myrs after the end of extension, these particles experience heating up to temperatures of 200–230 °C (Fig. 3C) as a result of overthrusting and syn-orogenic sedimentary burial (Fig. 1). Particle 8 is lying at middle-crust conditions (15 km, ∼ 380 °C) before onset of extension, and it represents the more distal part of the hyperextended domains considered in this study (Fig. 1). About 10 Myrs after initiation of rifting in Model A, the particle is heated to temperatures of 500–600 °C (Fig. 3C), corresponding to a local thermal gradient of 100–120 °C/km, as a result of rising of isotherms due to asthenospheric upwelling and drastic thinning of the lower crust. As extension increases and migrates toward the centre of the rift, upper crust and mantle become coupled. This leads to extensional unroofing and cooling of the particle in the footwall of a major detachment rooting into the mantle (Fig. 1). Finally, particle 9 shows a similar, albeit delayed by about 10 Myrs, synrift heating-cooling history (Fig. 3C). Overall, the rifting thermal history of particles 8–9, and most of the particles from the upper plate, is entirely preserved as the particles are not exposed to heating related to the subsequent nappe stacking evolution (Figs. 1 and 3C).

Predicted ages from Model A therefore reveal that the ZFT system is mostly insensitive to early orogenic processes. Only particles at crustal levels that record substantial thinning (middle crust) may provide insights on pre-collisional processes through partial resetting (path 8). In practice, this piece of information may be difficult to distinguish from the effect of heating during collision (paths 5 and 6), based purely on ZFT data. The two proximal margins do not record significant heating during thinning; they even cool down to temperatures below ZFT closure early during rifting. ZFT data from these parts of the margins are thus unable to resolve any pre-orogenic processes.

Model ZHe, AFT and AHe age distributions: syn-collisional cooling history?

Considering the very large ranges of values we modelled for sources of age variation, Figure 3 shows drastically different model ZHe, AFT and AHe signatures for particles of the upper plate and their counterparts in the lower plate, indicating distinct histories for the two plates. Patterns of thermochronological age predictions (ZHe, AFT and AHe) across the orogen are also very different from the model ZFT results.

Particles in the upper plate (paths 10–14) yield model ZHe, AFT and AHe ages that are generally pre-collisional (> 60 Ma), along with relatively long apatite MTL (lengths of > 13 μm account for 80% of the data set) and low predicted ZHe and AFT age variability (0–18%; Fig. 3B and Tab. SI03). These thermochronometric modelling results are broadly consistent with a scenario that involves (1) rapid cooling of the upper plate below the closure temperature of the AHe system during rifting and (2) limited to no heating of the plate above this closure temperature during collision. However, model AHe and AFT age patterns reveal some younging that is more pronounced for the hyperextended and the proximal domains (paths 10 and 14, respectively). Extremely younger AHe age predictions (0–19 Ma) and shorter model apatite MTL (the three shortest values for the upper plate, in the range 11.6–12.7 μm; Tab. SI03) for particle 10 reveal resetting by increasing sedimentary burial in the retroforeland basin in the last 20 Myrs (Figs. 1 and 3C). In contrast, the partial (AFT) to total resetting (AHe) of the model rifting information in apatite systems for particle 14 reflects two distinct phases of cooling (Figs. 1 and 3C): the first caused by tectonic denudation during extension (monotonic cooling through the apatite fission-track partial annealing zone over 25 Myrs; ca. 2 °C/Myr) and the second due to slower cooling thermal relaxation since the end of rifting (monotonic cooling through the apatite helium partial retention zone over 50 Myrs; ca. 1 °C/Myr).

Regarding particles in the lower plate (paths 1–8) and in the exhumed mantle domain (path 9), complexity manifested in model ZHe, AFT and AHe age patterns requires detailed analysis. Overall, ZHe, AFT and AHe ages modelled for upper crustal particles in the lower plate are younger than the ages predicted for particles with equivalent locations in the upper plate (Fig. 3A). Particles 5 (lower plate necking domain), 6–8 (lower plate hyperextended domain) and 9 yield ZHe, AFT and AHe age predictions younger than the onset of convergence (as old as 56 ± 3 Ma, 52 ± 3 Ma and 36 ± 6 Ma, respectively), along with long apatite MTL (13.5–14.6 μm; Figs. 3A and 3B). For most of these particles, model ZHe and AFT age ranges are separated by 5–10 Myrs, similarly to predicted AFT and AHe age ranges (5–15 Myrs, decreasing toward the lower plate), and model single-system data sets show a limited degree of age variability (8–40%, mostly 10–20%). Taken together, these thermochronometric modelling results seem to indicate a sustained cooling history of particles in the distal domains with relatively rapid crossing of both the higher-temperature ZHe and AFT closure windows due to convergence. Additionally, whatever the thermochronometric system, model single-system age patterns show fairly consistent younging from 50 to 10 Ma toward the lower plate, suggesting cooling migration from distal to proximal areas of the lower plate. As shown in Figure 1, these patterns reflect a sequence of accretion of lower plate crustal units and nappe stacking in the orogenic wedge. Close inspection of the time-temperature paths (Fig. 3C) and sequential evolution of the margins (Fig. 1) allows distinction between two different processes. Particles 5–7 record heating around 40–30 Ma (Fig. 3C), synchronously with syn-convergence tectonic and sedimentary burial (Fig. 1). This is followed by a period of cooling, starting when the weak middle-lower crust décollement is activated and thrusting and exhumation from mid-crustal levels initiate. In contrast, an earlier cooling history is observed for particles 8–9 (Fig. 3C), starting during rifting (75 and 65 Ma, respectively) and overall persisting at a sustained pace during the first 15 Myrs of convergence. Figure 1 highlights that cooling reflects a transient episode of coupling between crust and mantle prior to initiation of mid-crustal decoupling in the less thinned margin. This transient stage is a direct consequence of the rifting history that leads to the removal of the weak middle and lower crustal levels in the most distal portion of the margin. We infer that in distal domains the use of ZHe thermochronology, ultimately combined with AFT analysis, can yield constraints on early inversion mechanisms. In particular, we interpret the consistent and joint younging toward the lower plate of ZHe, AFT and AHe age predictions for particles in distal domains as intimately linked with the propagation of deformation from retro- to prowedge during convergence.

On the contrary, model ZHe and AFT age data sets for particles from the lower plate slightly stretched domain (paths 1–3) and its transition with the lower plate necking domain (path 4) show younging toward the orogen (Fig. 3A). They also present the highest degree of age variability obtained for the ZHe and AFT systems in this study (15–218%), spanning a time interval of 70–80 Myrs that encompasses both the rifting and convergence phases. A careful analysis of the time-temperature paths of particles 1–4 (Fig. 3C) reveals that the particles spent a significant amount of time in the ZHe partial retention (PRZ) and AFT partial annealing zones (PAZ) during and after rifting. In particular, each of the particles records heating in the last 40 Myrs. Timing and amount of heating are variable from one particle to another, as are timing and amount of sedimentary burial in the proforeland before final exhumation (Fig. 1), which determines whether the ZHe and AFT systems are fully or partially reset. For instance, particle 4 experiences heating up to 170 °C around 22 Ma (Fig. 3C), which places it in the ZHe closure window, far above AFT closure. Strikingly, the particle yields the second larger spread in ZHe age predictions in this study (19–82 Ma), with a narrow range of model AFT ages (15–18 Ma) and relatively long apatite MTL (13.1–13.6 μm) (Figs. 3A and 3B). Lastly, AHe ages predicted for the four particles are consistently much younger, although slight age variations illustrate different depths and temperatures reached by the upper crustal basement particles at the end of the model (Fig. 1).

Comparison with thermochronological age distributions in non-inverted, thermally-relaxed rifted margins (Model B)

As shown in Figure 4, model ZFT age distributions for non-inverted rifted margins (Model B) present no major difference with distributions obtained for the model of inverted rifted margins (Model A; Fig. 1). Most particles, especially those from the distal margins, indeed cool below 200 °C (path 9), even 140 °C (path 8), by the end of rifting in the two models (Figs. 3C and 4C). In details, we notice a reduced dispersion of model ZFT ages that reveals a simpler, continuous cooling history. Based on thermochronometric modelling results, this also manifests through overall longer predicted zircon MTL in Model B (Fig. 4B).

In contrast, model ZHe age patterns for Model B differ significantly from those for Model A, especially in the lower plate necking and hyperextended domains (paths 5–7) where pre-orogenic ages of ca. 80 Ma and older are remarkably dominant (Fig. 4A). In more distal domains, where particles experience initially greater temperatures (paths 8 and 9), postrift thermal relaxation results in cooling of the particles down to temperature below ZHe closure (140 °C) 5 to 15 Myrs after the end of extension (compared to 10 to 5 Myrs in Model A; Figs. 3C and 4C). Postrift cooling modelled in Model B is therefore overall slower. Yet, ZHe age predictions yielded by particles 8 and 9 in Model B (60–50 Ma; Fig. 4A) are similar in value to the true, early convergence ages obtained in Model A (Fig. 3A). Caution must therefore be taken when interpreting ZHe ages synchronous with early phases of convergence following rifting.

The effect of thermal relaxation identified earlier in model ZHe ages for particles 8 and 9 is even more pronounced in AFT age predictions, leading to slower crossing of the apatite fission-track PAZ by the particles (Figs. 3C and 4C) that manifests by much younger (34–20 Ma compared to 47–52 Ma) and more dispersed model AFT ages (71–149% compared to 11–10%; Figs. 3A and 4A), as well as much shorter predicted apatite MTL (10–13 μm compared to 13–14 μm; Figs. 3B and 4B). Such a thermal effect is also observed in the lower plate necking domain where the cooling particle 5 enters the AFT closure window around 5 Myrs after the end of extension in Model B (Fig. 4C), i.e. 35–40 Myrs earlier than it does in Model A (Fig. 3C). We therefore infer that interpreting independently ZHe and AFT age data sets as standalone evidence for an early orogenic event can be fraught with pitfalls.

AHe ages for Model B are all younger than predicted AFT ages, which is consistent with record of persisting regional postrift cooling by lower-temperature thermochronometers. Similar to AFT age predictions, model AHe ages yielded by particles 8, 9 and 5 in Model B are systematically younger than those predicted in Model A (Figs. 3A and 4A), reaching values of, or close to, zero. In the absence of syn-orogenic exhumation (Model A; Fig. 1), these particles indeed remain at temperatures above AHe closure (Fig. 4C). Nevertheless, the situation is different for particles originally lying at shallower upper crustal levels in the proximal and necking domains (3, 6, 7, 11, 12, 13 and 14) as they cool below 40 °C by the end of extension and overall yield consistent, pre-orogenic ages. In these domains, inversion and orogenic processes have no impact onto model AHe signatures, which offers the lowest-temperature thermochronometer considered in this study (AHe) great potential in confidently differentiating pre-convergence from syn-convergence histories.

Discussion

Key elements of syn-collisional cooling histories reflecting inherited rift architecture

Figure 5 synthesizes our main results. While the higher-temperature thermochronometer considered in this study (ZFT) documents the thermal history of particles during rifting, ZFT ages may reflect either crustal thinning in the lower plate necking and hyperextended domains (paths 5, 6 and 7) or postrift thermal relaxation in the more distal areas of the rifted margins (paths 8 and 9). In these latter distal areas, we also find that ZHe age predictions close in value to the timing of onset of convergence can reveal postrift cooling instead of collision initiation. Nevertheless, when model ZHe data are combined with AFT data that display strong signature contrasts between the suture zone and both margins, especially younger, more dispersed ages and shorter MTL (< 13 μm) in distal areas, we are able to diagnose unambiguously a control of postrift thermal relaxation onto the cooling history of particles. Close-to-zero model AHe ages in these regions are additional indicators of pure postrift thermal relaxation.

In contrast, near-identical model ZHe and AFT ages synchronous with early convergence in distal domains, along with narrow age variability in both data sets and long apatite MTL (> 13 μm), are indicative of tectonic inversion of the thin distal rifted margins. The older these ZHe and AFT age predictions, the closer the low-temperature thermochronology record to the timing of transition from extension to convergence and the location of particles to the upper plate. Although thermal signatures of particles in the lower plate necking and hyperextended domains are shown to correspond to distinct processes, predicted single-system age patterns skewing toward the lower plate, as displayed by both model ZHe and AFT data sets for our model of inverted rifted margins (Fig. 3A), provide insightful glimpses about timing of propagation of deformation from retro- to prowedge during convergence, including timing of early inversion mechanisms.

Application to the Pyrenees

The Pyrenean orogenic system, consisting of a prowedge on the southern, Iberian (lower) plate and a retrowedge on the northern, European (upper) plate (Fig. 6), is an exemplar setting to address crustal thermal and dynamic behaviour during orogenesis, as well as the role of rift inheritance into early orogenic processes using low-temperature thermochronology (Mouthereau et al., 2014; Vacherat et al., 2016; Ternois et al., 2019a, 2019b).

The Pyrenean orogeny (∼ 84–20 Ma) inverted a complex E-W- to ESE-WNW-trending Mesozoic rift system that formed in response to the opening of the Bay of Biscay associated with the migration northwards of the southern North Atlantic (Tugend et al., 2015; Barnett-Moore et al., 2016; Nirrengarten et al., 2017; Tavani et al., 2018; Angrand et al., 2020; Angrand and Mouthereau, 2021). Rifting developed mainly from late Aptian to early-to-middle Cenomanian (114–97 Ma) through extreme crustal thinning, associated with magmatic activity as well as heat advection and high-temperature-low-pressure metamorphism due to mantle exhumation (Ravier, 1957; Azambre and Ravier, 1978; Ravier and Thiébaut, 1982; Debroas, 1987, 1990; Golberg and Leyreloup, 1990; Lagabrielle and Bodinier, 2008; Lagabrielle et al., 2010; Clerc and Lagabrielle, 2014; de Saint Blanquat et al., 2016). Postrift lithospheric cooling and thermal subsidence from middle Cenomanian to latest Santonian was interrupted by the onset of N-S Iberia–Europe convergence, dated by the end of the Cretaceous Superchron at 83–84 Ma (Roest and Srivastava, 1991; Olivet, 1996; Rosenbaum et al., 2002a, 2002b; Schettino and Scotese, 2002; Macchiavelli et al., 2017). The major thermal pulse generated during rifting has been shown not to reequilibrate over the whole Pyrenean realm before the onset of convergence (Vacherat et al., 2014, 2016; Hart et al., 2017; Angrand et al., 2018).

Over the last 30 years, extensive low-temperature thermochronology studies in the Pyrenees have produced a world-class data set that shows diachronous development of the orogen, from east to west and from north to south, between the Late Cretaceous and the Miocene. A clear southward-younging thermochronological pattern has been emphasized, reflecting the prominent in-sequence tectonic growth of the Pyrenean orogen by accretion of Variscan basement units from the lower Iberian plate (e.g.Fitzgerald et al., 1999; Sinclair et al., 2005; Mouthereau et al., 2014; Fillon et al., 2020; Waldner et al., 2021). In contrast, the short delay (5–15 Myrs) between the end of rifting with associated high-temperature metamorphism and the onset of convergence makes it difficult to recognise the thermal signature of early orogenesis and to constrain early orogenic deformation (Vacherat et al., 2016; Ternois et al., 2019a, 2019b). Maximum temperatures recorded by numerical particles taken in the hyperextended and exhumed mantle domains are in line with those recorded in pre-orogenic Aptian-Cenomanian strata deposited over the inverted distal European margin (estimations of 350–500 °C, possibly 600 °C, obtained by Raman spectroscopy on carbonaceous material in the Metamorphic Internal Zone (Fig. 6); Vauchez et al., 2013; Clerc et al., 2015; Chelalou et al., 2016; Ducoux et al., 2019, 2021). Our integrated modelling approach can therefore provide a way to predict and help reconstructing the pre-orogenic rifted margin architecture and evolution in the Pyrenees based on low-temperature thermochronology signatures.

Figure 7 shows the distribution of ages documented across the Pyrenees by the four low-temperature thermochronology systems considered in this study. These ages have mainly been obtained from Variscan crystalline massifs composing the main thrust units of the Central Pyrenees (Fig. 6).

The North Pyrenean Zone is the retrowedge narrow north-verging fold-and-thrust belt comprising the inverted zones of Early Cretaceous extreme crustal thinning and mantle exhumation of the upper plate (Metamorphic Internal Zone; Fig. 6). In this domain of the Pyrenees, the central Arize and Trois Seigneurs basement massifs have yielded prerift to syn-orogenic bedrock ZFT and ZHe ages (169–71 Ma and 101–20 Ma, respectively), along with syn-orogenic bedrock AFT and AHe ages (50–35 Ma and 49–23 Ma, respectively) (Figs. 6 and 7; Vacherat et al., 2016). Fairly similar data are published from basement massifs and synrift strata of the western (Labourd-Ursuya − ZFT: 82 Ma; ZHe: 86–51 Ma; AHe: 49–35 Ma − Hart et al., 2017 and Mauleon − ZFT: 236–134 Ma; ZHe: 139–34 Ma − Vacherat et al., 2014), and eastern (Agly-Salvezines − ZFT: 75 Ma; ZHe: 117–57 Ma; AFT: 56–39 Ma; AHe: 55–40 Ma − Yelland, 1991; Gunnell et al., 2009; Ternois et al., 2019a) North Pyrenean Zone (Fig. 6). Altogether, these data are overall in line with age predictions for the upper plate exhumed mantle and hyperextended domains, where particles (paths 9 and 10) experience low to moderate crustal heating (up to 180–340 °C) during rifting then limited to no syn-orogenic burial (< 80 °C) (Model A; Figs. 1, 3 and 5). In this context, the model particles enter the zircon helium PRZ (particle 10), possibly the zircon fission-track PAZ (particle 9) during rifting, but may be not heated enough, or for long enough, for the ZFT or ZHe systems to be fully reset (e.g., crystals with high retentivity). This may thereby lead to preservation of prerift information by either or both systems (no reset), as well as record then preservation of synrift signatures (partial to total reset). Additionally, it is important to notice that particle 9, albeit cooling, experiences temperatures as high as 240 °C at the end of rifting, placing it within both ZFT and ZHe closure windows at the onset of convergence. Based on the argumentation in Ternois et al. (2019a), bedrock zircons with young, syn-orogenic ZFT or ZHe ages from the North Pyrenean Zone could be crystals with very low fission-track or helium retentivity. In this case, the crystal temperature sensitivity of these zircons could fall down to 180–200 °C (ZFT; Garver et al., 2005; Mesalles et al., 2014) or < 140–160 °C (ZHe; Guenthner et al., 2013), respectively, allowing at some point between rifting and convergence both reset of pre-orogenic information and record then preservation of syn-orogenic cooling (see sparse, sometimes debated, early convergence bedrock ZFT data by Yelland (1991) and ZHe data by Ternois et al. (2019a)). Overall, we infer that the thermochronological signature from the North Pyrenean massifs is indicative of the thermal evolution of the upper plate distal margin. This inference, based on our 2D thermally-coupled and dynamical modelling of low-temperature thermochronology data, is consistent with recently published pre-orogenic crustal architecture reconstructions for the Pyrenees (Mouthereau et al., 2014; Clerc et al., 2016; Ford et al., 2016; Teixell et al., 2016, 2018; Vacherat et al., 2016; Ternois et al., 2019a; Ford and Vergés, 2021; Waldner et al., 2021).

In the core of the Pyrenean range that consists of south-verging thrust sheets of Variscan metasediments and basement (the Axial Zone), all bedrock AFT and AHe age data from the northern Bethmale, Riberot and Marimana massifs are syn-orogenic (55–20 Ma) (Figs. 6 and 7). In contrast, both bedrock ZFT and ZHe data sets comprise a mixture of pre- and syn-orogenic ages (40–50 to 110–120 Ma). As such, these natural low-temperature thermochronology data fit well with our age predictions for the necking and hyperextended domains of the lower plate in Model A (crustal thickness of less than 7–5 km at the end of rifting), where particles (paths 7–8) record intense crustal heating (up to 600 °C) but limited syn-orogenic burial (< 200 °C) (Figs. 1, 3 and 5).

Both the Maladetta and Arties massifs, further south, have yielded ZFT and ZHe ages ranging 50–30 Ma and younger AFT and AHe ages (30–20 Ma) (Figs. 6 and 7). While younger ZHe ages are well reproduced in our model A for particles at the transition between the lower plate necking and hyperextended domains (paths 5 and 6; Figs. 1, 3 and 5), our modelling does not predict the young, syn-orogenic ZFT ages observed in the two Pyrenean basement massifs. It must be reminded that particles 5 and 6 are heated to maximum temperatures of 200–240 °C (∼ 9.3 km) due to syn-orogenic sedimentary and tectonic burial (Model A; Figs. 1, 3 and 5). This means that the particles enter the zircon fission-track PAZ but may not be heated enough, or for long enough, for the ZFT system to be fully reset (preservation of pre-orogenic ZFT ages). A slightly more important burial to 10 km (increasing maximum temperatures up to 300 °C) during orogenesis would permit total annealing of zircons, reset of pre-orogenic ZFT ages and record of syn-orogenic cooling by the ZFT system as is the case in the Pyrenean massifs. Alternatively, as argued earlier in this section, the zircons analysed from the two massifs could be crystals with very low fission-track retentivity, which lowers crystal temperature sensitivity and allows resetting of pre-orogenic information at lower temperatures. All in all, we conclude that observed ages from the Maladetta and Arties massifs most likely reflect cooling of the necking domain (up to 20-km-thick; Fig. 1) after a temperature increase during collision.

Finally, the southern Barruera and Bono massifs represent a domain of decreasing AHe ages (15 to 10 Ma; post-orogenic) and joint increase of AFT (20 to 40 Ma; syn-orogenic) and ZFT (from 50 to 60–120 Ma; syn- to pre-orogenic) ages toward the south, with scattered ZHe ages (50–20 Ma; syn-orogenic) (Figs. 6 and 7). These single-system age patterns are overall consistent with results obtained for particle 3 in our Model A, showing a generally sharp increase of AFT and ZHe ages toward the proforeland, along with large age variability in both single-system data sets (Figs. 1, 3 and 5). Based on our modelling, this suggests a thermal history involving denudation-related cooling of shallower crust during extension followed by exhumation then burial below syn-orogenic sediments deposited in the foreland basin (maximum temperatures of 140–170 °C). As shown in Figure 3C, old ZFT ages would most likely reflect temperatures inherited from the end of rifting, rather than burial and exhumation-related cooling from 300 °C during collision. We therefore infer that the thermochronological signature from the two Pyrenean massifs is indicative of the thermal evolution of the lower plate proximal margin.

Taken together, the above inferences not only imply that the upper plate distal domain has more chance to outcrop in the Pyrenees, north of the Axial Zone, but they also rule out any opportunity to reset appreciably the low-temperature thermochronology signal of the onset and early phase of orogenesis by burial. This therefore strengthens interpretations of Ternois et al. (2019a) about record of early inversion of the upper plate distal rifted margin by low-temperature thermochronology (mostly ZHe, AFT and AHe data) in the eastern retrowedge North Pyrenean Zone.

Conclusion

In this study, we perform detailed forward modelling to predict age patterns across a model tectonically-inverted hyperextended rift system for four widely used thermochronologic systems in orogenic settings (apatite (U-Th)/He, apatite fission-track, zircon (U-Th)/He and zircon fission-track systems). Our numerical low-temperature thermochronology data modelling shows that a combination of zircon (U-Th)/He (ZHe) and apatite fission-track (AFT) data allows diagnostic investigation of model orogen tectonics and offers the most valuable source of thermochronological information for the reconstruction of the crustal architecture of the model inverted rifted margins. Altogether, ZHe and AFT age data document the time-temperature history of rocks in a relatively large temperature range (40–250 °C), thus providing a record from which to unravel orogenic processes, variably spaced in time. Fission-track length distributions and intrasample age variability are as many additional constraints on the thermal history of a sample. Our modelling also demonstrates the need for caution when interpreting thermochronological ages as the evidence of exhumation.

Overall, the upper plate margin is predicted to preserve a clear rifting-related cooling signal in its proximal domains and early convergence cooling in its distal part. Postrift thermal relaxation is shown not to be sufficient for both the ZHe and AFT clocks to start during early convergence, clarifying the contribution of tectonic inversion to early convergence cooling of the upper plate distal margin. In contrast, thermal overprint due to later collision is expected within the lower plate margin. In this case, only ZFT data can provide constraints on the pre-collisional history but these will immutably document thermal histories during rifting.

We are aware that orogenic systems may involve complexities, in part inherited. These include, for instance, superficial decoupling levels such as evaporites that could delay the topographic response to convergence (Jourdon et al., 2020). However, the first-order similarities between our thermochronometric modelling results and low-temperature thermochronology constraints from the exhumed basement of the Pyrenees suggest that the initial geometry and thermo-mechanical properties of the basement of the margin are the main controlling parameters. In Alpine-type systems that involve a doubly-vergent orogen and result from the inversion of hyperextended rifted margins, we believe that our study can aid with the identification, especially during early orogenesis, of the different domains of a rifted margin that were accreted in the mountain belt.

Acknowledgments

This work was funded by the OROGEN research project (Total-BRGM-CNRS INSU). This manuscript greatly benefited from scientific discussion with the OROGEN team. We thank an anonymous reviewer and BSGF – Earth Science Bulletin Guest Editor O. Lacombe and Associate Editor P. Yamato for their constructive feedback and remarks. The authors, in particular ST, express gratitude to K. Gallagher for assistance with modelling, helpful insights and efforts put purposely into daily improvements of the software program QTQt.

Cite this article as: Ternois S, Mouthereau F, Jourdon A. 2021. Decoding low-temperature thermochronology signals in mountain belts: modelling the role of rift thermal imprint into continental collision, BSGF - Earth Sciences Bulletin 192: 38.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.