© C. Aumar et al., Published by EDP Sciences 2022

Highlights

• This study integrates in situ geochronological results with petrological and structural investigations to gain insight into the geodynamic evolution of deformed metamorphic rocks in the Agly Massif.

• Th-U/Pb dating of monazites and ⁴⁰Ar/³⁹Ar dating of muscovites and biotites yield a 94–127 Ma age bracket for the mylonitic deformation in the Agly Massif (Pyrenees, France).

• Based on the strain analysis in the Agly Massif, a purely divergent WNW-ESE trending continental rift significantly extended the lithosphere at the Iberia-European plate boundary during the Cretaceous.

1 Introduction

The geological province of the Pyrenees is part of the Variscan chain, which runs through the entire western European plate and was formed from the Devonian to the Carboniferous. Cadomian basement and Paleozoic cover have been strongly deformed and metamorphosed during the building of the chain.

Deeply eroded during the Permian, the Pyrenean province experienced a long period of calm in the Mesozoic, during which thick (up to 3000/4000 m) platform carbonates were deposited in conformity over the Triassic sedimentary sequence (e.g., Canérot, 1991), until an episode of extension in the Cretaceous at about 110–90 Ma. This tectonic event thinned the boundary between the Iberian and the European plates. Although coeval with a well-known HT-LP metamorphism of up to 600 °C in the Mesozoic sediments (e.g., Casteras, 1933; Ravier, 1959; Mattauer, 1968; Albarède and Michard, 1978; Golberg et al., 1986; Golberg and Leyreloup, 1990, Chelalou et al., 2016), this extension has long been considered as a minor episode associated with slight crustal thinning. Extensional structures have been described as a series of small pull-apart basins (Choukroune and Mattauer, 1978) resulting from the eastward rotation and/or translation of the Iberian Plate in the Cretaceous (e.g., Choukroune et al., 1973; Olivet, 1996; Sibuet, 2004; Tucholke et al., 2007; Jammes et al., 2009; Bronner et al., 2011).

Following the Cretaceous extension, the tectonic history of the Pyrenees area was concluded with the formation of the Pyrenees mountain chain (Fig. 1A). Due to the NNW movement of the Iberian microplate in the late Cretaceous and the Eocene (e.g., Choukroune, 1974), the Cretaceous basins were definitively closed and their tectonic inversion formed the Pyrenean chain (Teixell et al., 2018). Oriented east-west, this mountain belt deformed the Mesozoic and Cenozoic sediments. It also remobilized the old Paleozoic and Cadomian basement, often tilting it vertically, and it was crossed by large (i.e., several hundred meters) but local shear zones (e.g., Cochelin et al., 2017).

Studies of small peridotite bodies outcropping in the core of a few Cretaceous basins led to a reappraisal of the Cretaceous extension. Authors proposed that intense thinning of the crust took place in the basins, thus explaining the exhumation of the underlying mantle (e.g., Lagabrielle and Bodinier, 2008; Jammes et al., 2009). Since then, careful studies have shown that the Cretaceous extension had probably been underestimated, and that this tectonic event involved strong ductile deformation and major thinning of the continental lithosphere (e.g., Clerc et al., 2012; Vauchez et al., 2013; Clerc and Lagabrielle, 2014; Lagabrielle et al., 2016, 2019a, 2019b; de Saint Blanquat et al., 2016; Corre et al., 2016; Asti et al., 2019).

According to these authors, the Cretaceous event can be described as an east-west trending continental rift (i.e., north-south extension). In this hypothesis, it is supposed that thinning of the crustal basement led in some places to its total disappearance (i.e., lateral extraction), forming necking zones, which would have put the underlying mantle into contact with the pre-extensional Mesozoic cover, explaining the HT-LP metamorphism in the latter. A boudin of crustal basement would have separated each necking zone.

In the rifting model, crustal necking zones underwent very high strain. Opinions diverge on the behavior of interleaved boudins. Depending on the authors, they are considered to be partially ductile, to undergo internal strain (e.g., Clerc and Lagabrielle, 2014), or to be brittle, behaving as rigid blocks (e.g., Duretz et al., 2019). In the first hypothesis, the basement would have deformed to create a pinch and swell structure. In such structures, the swells are simply less deformed than the pinches and they undergo internal strain.

It should be noted that previous structural and geochronological studies suggest that part of the ductile deformation recorded in the Paleozoic basement of the North Pyrenean Zone might have resulted from a Cretaceous event dated between 110 Ma and 90 Ma (Passchier, 1984; de Saint Blanquat et al., 1986; Costa and Maluski, 1988; Monié et al., 1994; Odlum and Stockli, 2020). If correct, this Cretaceous deformation of the basement would be coeval with the HT/LP metamorphism recorded in the pre-rift Mesozoic sediments (see review in Clerc et al., 2015).

To check the pinch and swell hypothesis, it seems pertinent to study a Paleozoic Massif considered in the rifting model as a tilted rigid block bounded by two necking zones (e.g., Duretz et al., 2019). The Agly Paleozoic Massif, in the eastern part of the North Pyrenean Zone, is well suited to such a study because the pre-Mesozoic basement is well exposed in the area and is thought by most authors to have undergone mainly Variscan deformation, metamorphism and magmatism (e.g., Tournaire Guille et al. 2019; Siron et al., 2020). In this study, the focus has been placed on the deeper strata of the massif, which may have been deformed during the Cretaceous rifting. Determining the age of the main deformation observed in this lower unit should provide a better understanding of the mechanics of the rifting process.

Fig. 1 (A) Simplified geological map of the North Pyrenean Zone (NPZ) and location of the Agly Massif in the Pyrenean belt. (B) Geological map (modified from Fonteilles et al., 1993 and Berger et al., 1993) showing the location of the Agly Variscan Massif between the Saint Paul de Fenouillet and Bas Agly sedimentary basins to the north and the Boucheville sedimentary basin to the south. The North Pyrenean Fault (NPF) separates the Northern Pyrenean Zone from the Axial Zone. Stretching lineations from Delay (1989), Bouhallier et al. (1991), Choukroune et al., (1990), Vauchez et al., (2013) and this study.

2 Geological setting

From south to north, the Pyrenean chain is subdivided into three zones: the South Pyrenean Zone, the Axial Zone and the North Pyrenean Zone. The Axial Zone is separated from the North Pyrenean Zone by a major fault called the North Pyrenean Fault (NPF). The NPF results from the re-activation during the Pyrenean orogeny of the transform fault, which separated the former Iberian microplate from Europe.

The Agly Massif (Eastern Pyrenees, France) is the easternmost Paleozoic Massif of the North Pyrenean Zone (Fig. 1A) and belongs to a series of massifs that outcrop along this zone. It is bounded by the Cretaceous Boucheville basin to the south and by the Cretaceous Agly and the St-Paul-de-Fenouillet basins to the north (Fig. 1B), which underwent tectonic inversion in the Eocene. These basins contain both pre- and syn-rift Mesozoic sediments. The pre-rift sediments were deposited directly onto the eroded Variscan basement and correspond to a classic carbonate platform including limestones, dolomites and marls. At their base, a thick level of Upper Triassic marls and evaporites (Keuper deposits) is generally observed, whose thickness varies from 1000 m to 2700 m (e.g., Jammes et al., 2010). Locally, up to 5000 m of Albo-Cenomanian black flysch were deposited during the Cretaceous extension, highlighting the importance of subsidence during this period.

These three sedimentary basins are considered by the proponents of the rifting model as being crustal necking zones (e.g., Vauchez et al., 2013; Clerc et al., 2016) where lateral extraction of the basement put the uppermost lithospheric mantle in contact (or nearly so) with overlying Mesozoic sediments, generating HT metamorphism of up to 600 °C in the sediments (e.g., Ravier, 1959; Golberg and Leyreloup, 1990). Peak temperature estimates in Mesozoic sediments of the Boucheville basin south of the Agly Massif are up to 580 °C at a depth of about 6/7 km (Chelalou et al., 2016). U-Pb geochronology and field observations reveal that magmatic sills intruded the sediments at 97 Ma (Chelalou et al. 2016). According to these authors, the basement below this sedimentary basin was probably very thin, or even absent, due to extreme thinning. To the north of the Agly Massif, the Mesozoic sedimentary series of the Agly basin recorded a thermal episode up to 390 °C at 100 Ma, associated with top-to-the-NE mylonitization (Vauchez et al. 2013). Thinning of the underlying basement in this northern basin is considered to have been less significant than in the south.

Generally described as a dome slightly inclined towards the east, the Agly Massif is not thought to have been affected by the late Pyrenean orogeny, thus preserving the Variscan deformation (Olivier et al., 2004). It can be divided into a Lower Gneiss Unit (LGU) and an Upper Schist Unit (USU) (Fig. 1B) (Fonteilles, 1970). Since the USU is Paleozoic in age, the age of the LGU protolith could be either Paleozoic or Cadomian, or both as indicated by recent U-Pb zircon ages (Tournaire Guille et al., 2019). The two units are separated from each other by the Latour-Caladroy shear zone, interpreted as a major top-to-the-NNE detachment fault (e.g., Delay, 1990; Bouhallier et al., 1991).

A HT-LP metamorphic episode is recorded in the Agly Massif and generally attributed to the late-orogenic Variscan extension, during the late Carboniferous (Fonteilles, 1970; Bouhallier et al., 1991; Olivier et al., 2004; Siron et al., 2020). Metamorphic parageneses point to PT conditions of 600/700 °C at a depth of about 10/12 km (Guitard and Raguin, 1958; Fonteilles, 1970; Vielzeuf, 1984; Olivier et al., 2004; Siron et al., 2020). In connection with this high geothermal gradient, the whole series of the Agly Massif is intruded by late Variscan granitic plutons (i.e., 310/300 Ma) (Fig. 1B) (Michard and Allegre, 1975; Respaut and Lancelot, 1983; Vielzeuf and Pin, 1991; Fourcade and Javoy, 1991; Delaperriere and de Saint Blanquat, 1994; Althoff et al., 1994; Olivier et al., 2008; Tournaire Guille et al., 2019), which indicates that partial melting of the crust occurred at depth at that time.

The LGU consists of paragneisses interleaved with lenses of orthogneiss. Compositions range from quartzo-feldspathic gneisses with biotites and muscovites to garnet-sillimanite-cordierite kinzigites. With the notable exception of a 4 km-wide EW corridor in the northern part of the Massif, the LGU displays flat-lying foliations, which have always been considered to be relics of the Variscan deformation. Gneisses are well deformed and often show mylonitic foliations with a uniform N10E/N30E stretching lineation (Delay and Paquet, 1989; Bouhallier et al., 1991; Vanardois et al., 2020). It is generally assumed that the post-Variscan deformation simply led to foliations being tipped into a vertical orientation during the Pyrenean orogeny.

The aim of this study is to determine the age of the main deformation observed in the LGU of the Agly Massif in order to clarify whether it resulted from the Variscan or the Cretaceous tectonic event. Given that the deformation in the Massif is very heterogeneous, we mainly collected samples in flat-lying highly deformed zones with a well-defined N20E oriented stretching lineation When strain gradients occurred over only a short distance, we also collected samples from less deformed zones. The Table 1 indicates the GPS coordinates and mineralogical description of each sample while the Figure 1B shows their location on the geological map.

3 The main pervasive deformation (MPD) in the LGU

Foliations are predominantly flat-lying and concordant with lithological units. This foliation is ubiquitous regardless of the strain intensity. The deformation is heterogeneous, with strain gradients on all scales generating mylonitic foliations (Fig. 2A). In these zones, strain can be strong enough to render minerals too small to be observed both in the field and under a microscope (i.e., ultra-mylonites).

Quartz ribbons, biotite clusters or elongated minerals like feldspar or garnet are oriented in a direction that defines the stretching lineation, considered as being the long axis of the strain ellipsoid (i.e., the λ₁ axis). As already noted by previous authors (Delay and Paquet, 1989; Bouhallier et al., 1991; Vanardois et al., 2020), it has an average orientation of N20E (Fig. 1B). This marked, omnipresent stretching lineation, associated with mylonitization, defines the main pervasive deformation (MPD) imposed on the LGU.

The MPD is associated with numerous shear bands, which are slightly oblique to the foliation plane. In mylonites and ultra-mylonites, shear bands and foliations are generally parallelized to form a coherent MPD environment. Shear senses can be determined using kinematic indicators like shear bands, the sigmoidal shapes of recrystallized tails around feldspars or garnets, asymmetric pressure shadows and feldspars with domino-like structures (e.g., Simpson and Schmid, 1983; Hanmer and Paschier, 1991) (Figs. 2B–2D).

Two opposite senses of shear can be observed in the same areas, both associated with the NNE-SSW stretching lineation. These two shear senses, top-to-the-SSW and top-to-the-NNE, were also noted by Bouhallier et al. (1991) and Vanardois et al. (2020). The question is whether these opposing senses of shear correspond to superimposed deformations or to a single deformation episode.

Despite some rare low-strain zones, it is difficult to find evidence of an earlier foliation in the field. This is due to the obliquity between the MPD foliation and the coeval numerous shear bands, which prevents any clear evidence of a previous surface. In high-strain zones, all surfaces are parallelized, including any possible earlier foliation. However, this earlier foliation can be seen in some thin sections where it is strongly deformed into isoclinal folds (Fig. 2E). Only the hinges are detectable while flanks are sheared and parallel to the new MPD foliation. This scarce relictual foliation (RF) shows that the strain intensity of the MPD was strong enough to erase almost all previous deformations.

Fig. 2 Pictures A, B, C and D refers to the Figure 1 where the location of the photos is shown (orange circle). (A) MPD flat-lying mylonitic shear zone in orthogneiss (east of Mas Chiffre-Rnes). (B, C and D) Kinematic indicators of the MPD (St-Martin/Cuxous area). (E) Thin-sections showing the folding of the relictual fabric (RF) transposed onto the MPD foliation.

3.1 Strain pattern of the MPD

Two cases can be described when looking at kinematic indicators in thin section. In the first one, most thin sections show a clear predominant sense of shear, either top-to-the-SSW or top-to-the-NNE. Both shear directions are strictly parallel to the NNE-SSW oriented stretching lineation. Shear bands are shallowly dipping (10 to 30 degrees to the SSW or the NNE) and clearly oblique to the horizontal foliation. Conjugate shear bands may be present but the predominance of one set over the other results in a structural asymmetry that indicates a non-coaxial deformation on the scale of the thin section. This confirms observations made in the field.

The second case corresponds to anastomosing shear bands, which isolate lenses whose long axis is parallel to the flat-lying foliation (Fig. 3). Each individual lens is limited along its two upper faces by a top-to-the-SSW sense of shear and a top-to-the-NNE sense of shear. Likewise, the two lower faces are limited by top-to-the-NNE and top-to-the-SSW shearing. This symmetric geometry is especially well expressed around feldspar porphyroclasts (Figs. 3B and 3C). At the scale of the whole thin section, asymmetry may dominate according to the situation with respect to the lens. Predominant top-to-the-SSW sense of shear is observed along two opposite faces of the lens while top-to-the-NNE sense of shear is dominant along the two others faces. For instance, in the Figure 3A, top-to-the-SSW sense of shear (SC structures) is well visible in the upper half of the thin section while top-to-the-NNE sense of shear (SC structures) prevails in the lower half. We interpret this thin section as the SSW ending of a larger lens. The shear bands surrounding the lenses exhibit mylonitic or ultra-mylonitic deformation whereas the core of the lenses is significantly less deformed. The formation of lenses results in a more or less symmetric structural pattern on the scale of the thin section (Fig. 3A).

It is known that the strain regime can be deduced from the degree of symmetry of the structural pattern (Choukroune et al., 1987). Structural symmetry involving conjugate sets of shear zones is the result of bulk inhomogeneous coaxial strain (e.g., Choukroune and Gapais, 1983; Gapais et al., 1987). This statement helps to interpret the lens-shaped pattern in thin section. According to the flat-lying foliation and the orientation of the long axis of lenses, this pattern accounts for coaxial vertical shortening.

This observation strongly suggests that top-to-the-SSW and top-to-the-NNE senses of shear in the Agly Massif do not correspond to superimposed deformations. Both sets of shear bands, which display the same orientation but opposite senses of movement, are likely to be coeval, resulting in the formation of the lenses. Also, the constant NNE-SSW direction of the stretching lineation and the two opposite shearing senses, as measured in the field, probably reflect a 2D deformation and are consistent with bulk coaxial plane strain.

Surprisingly, lenses are not commonly seen at outcrop scale (i.e., meter-scale). According to the Von Mises criterion, ductile conjugate shear zones initiate at 45° to the flattening direction. With increasing strain, they rotate towards the flattening plane, which increases the angle of the shortening direction (i.e., ≈ 120°).

A plausible explanation of the lack of macroscopic lenses could be that the relictual foliation acted as a strong prior mechanical anisotropy, preventing the formation of conjugate shear planes intersecting it at a high angle. The energy needed to initiate conjugate shear zones at 45° to the previous well-developed flat-lying foliation plane (i.e., the relictual foliation) during bulk vertical shortening is probably too great compared to that needed for opposite shear zones at a low angle to the foliation plane.

In thin-section, microscopic lenses were initiated on resistant heterogeneities like feldspar and garnet porphyroclasts (Figs. 3B and 3C). Although rare, these small lenses can sometimes be observed in the field on a centimeter/decimeter scale in augen orthogneisses, probably due to the presence of large feldspars and garnets (Fig. 3E). Such heterogeneities do not exist on a larger scale in the LGU. Based on this observation, the deformation in the LGU took place through the superposition of slightly oblique shear zones with opposite senses of shear, which flattened over time due to bulk coaxial shortening to become more or less parallel. This unusual, staged geometry can be observed in a single thin section (Fig. 3E).

Following the MPD, brittle deformation is also expressed in thin section by conjugate Coulomb-type fractures dipping at 60 degrees (Figs. 3D, 4E and 4F). They either crosscut the whole thin-section, or deformational structures like quartz ribbons and shear bands, or even individual minerals like feldspar and garnet. According to the Mohr–Coulomb theory, the maximum principal stress of the stress ellipsoid (i.e., the σ₁ axis) is the bisector of the 60° acute angle. The dip of these conjugate brittle faults then accounts for vertical shortening. Likewise, given the azimuth of these fractures, the minimum compressive stress (σ₃) is horizontal and oriented NNE-SSW.

Fig. 3 Thin sections showing (A) dominant top-to-the-SSW senses of shear (SC structures) at the upper half of the thin section and dominant top-to-the-NNE senses of shear at the lower half interpreted as the SSW ending of a microscopic lens, (B and C) symmetric MPD lenses around resistant feldspars, (D) detail of a late semi-brittle normal fracture crosscutting the whole thin section, (E) MPD domino-type bookshelf feldspars indicating parallel but opposite senses of shear. The thin section consists of parallel bands indicating alternating top-NNE and top-SSW shear senses. In each parallel band, bookshelf feldspars are consistent with the flow revealed by quartz fabrics, pressure shadows and shear bands, (F) outcrop showing centimetre-scale MPD lenses (Saint Martin area). Picture E refers to the Figure 1 where the location of the photos is shown (orange circle).

Fig. 4 Photomicrographs of quartzofeldspathic gneisses with biotites and garnets. (A) Almond-shaped K-feldspar surrounded by fine-grained aggregates of biotites and recrystallized quartz ribbons. (B) Almond-shaped K-feldspar later deformed into a domino-type. (C) Domino-type feldspar later partially ground at its tails. (D) Almond-shaped K-feldspar later deformed into a domino-type and finally completely disrupted in a cataclastic manner at its tails. The feldspar is surrounded by fine-grained quartz ribbons. (E) K-feldspar crosscut by conjugate Coulomb-type fractures dipping at 60 degrees. (F) Fractured garnet crosscut by conjugate Coulomb-type fractures dipping at 60 degrees.

3.2 Temperature conditions of the MPD

In thin sections of the least deformed samples, the primary metamorphic assemblage including feldspar, garnet, sillimanite, cordierite, biotite and muscovite is still recognizable (Fig. 2C), but is partially replaced by a syn-kinematic assemblage including biotite and muscovite.

In highly sheared gneisses, which display a very strong NNE-SSW stretching lineation, the primary metamorphic or magmatic assemblage is fully replaced by a new assemblage made of recrystallized fine-grained quartz and feldspar and syn-kinematic micas. When present, garnet appears as undeformed crystals or is sheared by microfracturing. Sigmoidal shapes resulting from garnet fracturation are consistent with the sense of shear giving an illusion of plasticity (Fig. 2D). In contrast, almond-shaped K-feldspar and plagioclase indicating plastic deformation are observed (Figs. 4A and 4B, samples AG17-J5b and AG17-J6b). The quartz has been dynamically recrystallized and forms aggregates of fine, equant grains, which in turn have undergone ductile deformation.

However, in the same thin section, feldspar can also behave in a “domino” fashion by activation of one of its two cleavages (Figs. 3E, 4B–4D). Together with dynamic quartz recrystallization, the formation of quartz ribbons is observed in narrow bands of very high strain (Figs. 3E, 4B–4D).

In thin sections it can be seen that the outer rim of feldspars is often fragmented into small angular pieces (Figs. 3C, 4C and 4D). The entire feldspar can also be ground up and fragments dispersed and re-oriented along the shear bands. Likewise, some of the quartz tails around feldspar crystals have been completely disrupted in a cataclastic manner forming a fine-grained, crushed matrix in many shear bands. As already described in the “strain pattern” section, the final stage of deformation corresponds to conjugate brittle faults cross-cutting the foliation plane as a whole and, on a small scale, feldspar and garnet crystals (Figs. 4E and 4F).

Although mineral behaviour also depends on pressure and water content, thin section analysis gives rough but coherent estimates of the thermal conditions at the time the shearing deformation, associated with a well-defined N20E stretching lineation, took place. The two dominant minerals in the LGU are quartz and feldspar. The way these two minerals have deformed in shear bands or pressure shadows provides valuable evidence about the temperature at the time of deformation.

Almond-shaped K-feldspars indicate ductile deformation of that mineral. Feldspar, because of its very low symmetry, is highly resistant to plastic flow. Plastic deformation of feldspar implies recovery, and recrystallization can only occur above 500 °C, thus in the amphibolite facies (e.g., White et al., 1975; Voll, 1976).

Domino-type feldspars occur below the field of plastic flow for that mineral, generally in the range of 300 °C to 500 °C (e.g., Gapais, 1989), thus at a temperature significantly lower than for the almond-shaped feldspars. Narrow bands of quartz ribbons are typical of strain localization. Together with domino-like feldspar structures, they indicate temperatures below amphibolite facies, at less than 500 °C. This is in agreement with the alteration of feldspar into syn-kinematic chlorite and sericite, as frequently observed in thin section.

Finally, tails of feldspars disrupted and grounded in a cataclastic manner along shear bands, as well as late conjugate brittle faults cross-cutting thin-sections or individual porphyroblats, reveal a stage of brittle deformation at a temperature below 300 °C.

The crucial point is that almond-shaped or domino-like feldspars, dynamic recrystallization of quartz and quartz ribbon formation, as well as cataclasis creating angular clasts in a fine-grained matrix, can all be observed associated with shear bands having the same shear sense in individual thin sections. Together with the late conjugate brittle faults, we conclude that all these structural features formed during the same continuum of deformation, which therefore occurred under decreasing temperature conditions. The MPD started with ductile strain at relatively high temperatures (500 °C) and ended with brittle deformation at a very low temperature (150/200 °C).

4 Monazite LA-ICPMS dating

This study focuses on three distinct areas of the LGU where samples have been collected for petro-structural and geochronological (Th-U/Pb and ⁴⁰Ar/³⁹Ar) analyses (Fig. 1B). Sample names, GPS coordinates, mineral and deformation characteristics as well as isotopic ages are reported in Table 1.

Th-U/Pb analyses of monazite have been performed by LA-ICPMS on thin sections. A laser spot size of 9 µm was used for ablation. Analytical methods are presented in detail in Supplementary Material S1. All Th-U/Pb data are quoted at the 2σ level and reported in Table S2 (Supplementary Material S2). Data are projected on ²⁰⁶Pb/²³⁸U versus ²⁰⁸Pb/²³²Th (Figs. 6A, 6C, 6E, 7A, 7C and 7E) diagrams as well as in Tera–Wasserburg (Figs. 6B, 6D, 6F, 7B, 7D and 7F) or Concordia diagrams (Fig. 6F), as commonly used to present in situ monazite analyses. This method provides the best results given the behaviour of the three isotopic systems used for monazite dating. Except for four concordant inherited ages between 500 Ma and 596 Ma (samples AG17-J6b and AG17-J15b), all concordant ²⁰⁸Pb/²³²Th vs. ²⁰⁶Pb/²³⁸U ages (and ²⁰⁷Pb/²³⁵U-²⁰⁶Pb/²³⁸U) are scattered between 364 Ma and 104 Ma with two broad peaks at ca. 280–320 Ma and 105–120 Ma. This bimodal distribution suggests the existence of two main generations of monazite, referred to below as M1 (Variscan) and M2 (Cretaceous). Intermediate ages are highly scattered and probably reflect mixing.

Fig. 5 BSE images of monazite (samples AG17-A2, AG17-J15b and AG17-J6b). Monazites are found as inclusions within large garnets (A) or found in the matrix (B, C and D) in contact with biotites (D), K-feldspars (B) or garnets (C). Large monazite grains are surrounded by small grains or recrystallized branches (C) or can be clearly deformed (B). Small monazite grains can also form clusters aligned along the foliation plane (D). La (left) and Y (right) X-ray maps of monazites (E) in matrix and (F) at the edge of garnet in sample AG17-J15a. Circles show the location of the LA-ICPMS pits (9 μm) and their corresponding 208Pb/232Th ages. Circles show the location of the LA-ICPMS pits (9 μm) and their corresponding 208Pb/232Th ages (In red: concordant Cretaceous age; in blue: mixing age; in black: Variscan age).

Fig. 6 (A) 208Pb/232Th versus 206Pb/238U diagram for monazites in sample AG17-A2. Dashed ellipses correspond to discordant data. Dark ellipses correspond to concordant ages affected (or not) by mixing between M1 and M2 monazite domains. White ellipses correspond to M2 analyses affected by common Pb contamination and mixing with M1 domain. Grey ellipses correspond to analyses used for the Discordia calculation. (B) Tera–Wasserburg diagram for the same sample. (C) 208Pb/232Th versus 206Pb/238U diagram for monazite in sample AG17-A3. (D) Tera–Wasserburg diagram for the same sample. (E) 208Pb/232Th versus 206Pb/238U diagram for monazite in sample AG17-H1. (F) Concordia diagram for the same sample.

Fig. 7 (A) 208Pb/232Th vs. 206Pb/238U diagram for monazites in sample AG17-J15a. Dashed ellipses correspond to discordant data. Dark ellipses correspond to concordant ages affected (or not) by mixing between M1 and M2 monazite domains. White ellipses correspond to M2 analyses affected by common Pb contamination and mixing with the M1 domain. Grey ellipses correspond to analyses used for the Discordia calculation. (B) Tera–Wasserburg diagram for the same sample. (C) 208Pb/232Th vs. 206Pb/238U diagram for monazites (sample AG17-J15b) and (D) associated Tera–Wasserburg diagram. (E) 208Pb/232Th versus 206Pb/238U diagram for monazites in sample AG17-J6b. Dashed ellipses correspond to discordant data. Grey ellipses correspond to analyses used for the Discordia calculation. This sample is characterized by two populations of chemically distinct monazites, which allow two Discordia lines to be drawn. See the text for detailed explanations. (F) Tera–Wasserburg diagram for the same sample.

4.1 Textural and chemical characteristics of the monazite grains

In all samples, which mainly range from proto- to ultramylonitic gneisses, monazite occurs both in the matrix and as inclusions (20 µm to 50 µm) primarily in garnet and occasionally in feldspar porphyroblasts (Fig. 5A). Matrix monazites are commonly subhedral-elongated grains (50 µm to 300 µm), with irregular rims and cracks (Fig. 5B). In ductile shear bands or mylonitic zones they show specific textures such as ramifications, small satellite grains (5 µm to 10 µm) surrounding larger monazite grains (20 µm to 100 µm) (Figs. 5B and 5C) or clusters of small grains (10 µm) elongated parallel to shear bands (Fig. 5D) and characteristic of dissolution-recrystallization processes. These monazites are commonly associated with secondary minerals such as allanite, apatite, xenotime and small biotite grains (5 µm to 10 µm). In garnet fractures, monazite displays recrystallization textures and/or an elongated shape along the fracture direction and is associated with chlorite, muscovite and biotite (Fig. 5C).

Variations in the chemical composition of monazite (Tab. 2) are dependent on the mineral assemblages and the textural position of the monazite grains, as well as their age. Two generations of monazite are observed, Variscan (M1) and Cretaceous (M2), whose chemical composition differs from garnet-bearing rocks to garnet-free rocks.

In garnet-bearing rocks, M1 monazites included in garnet without destabilization textures have a homogeneous and low Y content (0.5–1.5 wt%, Tab. 2). Most of the M1 matrix monazites display patchy Y zoning, with Y-rich and Y-poor domains. M2 monazites either have a thin (1–10 µm) Y-rich (3–4 wt%) rim around an Y-patchy M1 core, or exist as satellite grains and clusters usually with a high Y content (3–3.5 wt%) (Figs. 5E and 5F). M1 monazites are Th-rich (4–7 wt%) compared to M2 monazites, which have a lower Th content (1–3 wt%). HREE (Gd + Sm) content is higher in M2 monazites than in M1, while LREE (La+Ce+Pr+Nd) content is the same (Figs. 5E and 5F and Tab. 2).

In garnet-free rocks, M1 monazites show a heterogenous Th and LREE distribution (Tab. 2). Large cores of M1 display patchy zoning characterized by variations in Th (4–7 wt%) and LREE (45–56 wt%) content and a high, homogeneous Y content (3–4 wt%). M2 monazites correspond to a Y-poor (0.5–1.5 wt%) thin rim around the large M1 core with a heterogeneous Th (0–5 wt%) and LREE (50–55 wt%) content. Some of the M2 monazites display a totally distinct chemical composition characterized by a very low Th (from undetectable to 1 wt%) and high REE (60–63 wt%) content (Tab. 2). The Y content (0.5–1.5 wt%) is the same in the Th-poor and Th-rich monazites.

4.2 Th-U/Pb geochronological results

Most of the Th-U/Pb analyses measured in M2 monazite domains (Cretaceous) are discordant due to various proportions of common Pb contamination and variable amounts of mixing with the M1 monazite domains. This results in analytical points plotting above the Concordia in the Tera–Wasserburg diagram (Figs. 6B, 6D, 7B, 7D and 7F), and below the Concordia in the U/Pb Concordia and the U-Th/Pb diagrams for most of the monazites analysed. In both diagrams, analytical points affected by common Pb contamination alone plot along a mixing line between the common Pb isotopic composition (reported in Supplementary Material S2 and consistent with the common Pb isotopic composition calculated by Stacey and Kramers, 1975) in the upper intercept, and the monazite age in the lower intercept. In sample AG17-J6b, some of the monazites are characterized by extremely low Th contents and more abundant common ²⁰⁸Pb than radiogenic ²⁰⁸Pb. In the Th-U/Pb diagram, analytical points of Th-poor monazites plot below the Concordia curve along a mixing line between the common ²⁰⁸Pb/²³²Th ratios in the “lower” intercept and the monazite age in the “upper” intercept (Fig. 7E).

In sample AG17-A2, concordant analyses (40 analyses from 6 monazite grains) in the ²⁰⁶Pb/²³⁸U vs. ²⁰⁸Pb/²³²Th diagram (Fig. 6A) spread from 351 Ma to 115 Ma. Analyses of the M1 monazites are mostly concordant, while analyses of the M2 monazites are usually discordant due to common Pb contamination and mixing. Six of the M2 monazite analytical points plot along a mixing line between the common Pb composition and the lower intercept age at around 109 Ma (MSWD = 1.8). One of these analyses yields a ²⁰⁶Pb/²³⁸U–²⁰⁸Pb/²³²Th concordant age of 108 ± 4 Ma. In the Tera–Wasserburg diagram (Fig. 6B), five concordant analyses spread from 349 Ma to 261 Ma. Ten discordant analyses are used to draw a Discordia line, which intercepts the Concordia curve at 111 ± 3 Ma (N = 10, MSWD = 2.0). In sample AG17-A3, 65 analyses from 15 monazite grains are plotted in the ²⁰⁶Pb/²³⁸U vs. ²⁰⁸Pb/²³²Th diagram (Fig. 6C) and in the Tera–Wasserburg diagram (Fig. 6D). In both diagrams, the younger the ages, the more discordant the analyses (Figs. 6C and 6D). In the ²⁰⁶Pb/²³⁸U vs. ²⁰⁸Pb/²³²Th diagram (Fig. 6C), concordant analyses are between 298 Ma and 157 Ma. It was not possible to define a Discordia line in this diagram. In the Tera–Wasserburg diagram, most of the analyses are discordant due to common Pb contamination and mixing. Concordant analyses are between 161 Ma and 287 Ma. Four of the analyses from the discordant analytical points in M2 monazites define a Discordia line with a lower intercept age at 112 ± 15 Ma (MSWD = 3.2), while the upper intercept is consistent with the common Pb isotopic composition calculated at the same age (common Pb isotopic composition at 110 Ma: ²⁰⁶Pb/²⁰⁴Pb = 18.536; ²⁰⁷Pb/²⁰⁴Pb = 15.621; ²⁰⁸Pb/²⁰⁴Pb = 38.43; Corresponding age: 4995 ± 250 Ma, Supplementary Material S2).

In sample AG17-H1, 30 measurements were performed on a set of 15 monazites and plotted in the ²⁰⁶Pb/²³⁸U vs. ²⁰⁸Pb/²³²Th diagram (Fig. 6E) and in the Concordia diagram (Fig. 6F). All except one data point are concordant in the Concordia diagram, yielding ages between 351 Ma and 291 Ma in both systems. No common Pb contamination is observed.

In sample AG17-J15a, 39 analyses were obtained from 8 monazite grains. In the ²⁰⁶Pb/²³⁸U vs. ²⁰⁸Pb/²³²Th diagram (Fig. 7A), ages are mostly concordant between 106 Ma and 321 Ma. In the Tera–Wasserburg diagram (Fig. 7B), concordant ages are between 103 Ma and 319 Ma. Three discordant analytical points plus the concordant data align at 103 Ma, all obtained from the M2 monazites, making it possible to define a Discordia line which intercepts the Concordia at 103 ± 9 Ma (MSWD = 2.8, upper intercept in Supplementary Material S2). In sample AG17-J15b, 42 analyses were provided by six monazite grains. In the ²⁰⁶Pb/²³⁸U vs. ²⁰⁸Pb/²³²Th diagram (Fig. 7C), the ages are mostly concordant between 104 Ma and 364 Ma. In the Tera–Wasserburg diagram (Fig. 7D), concordant analyses are between 122 Ma and 282 Ma, with one older domain at 507 Ma in a matrix grain. The younger the ages, the more discordant the data. Four of the discordant analytical points measured in the M2 monazites plus the concordant data point at 104 Ma plot along a mixing line between the common Pb composition and the lower intercept age at 102 ± 6 Ma (MSWD = 2.8, upper intercept in Supplementary Material S2).

Eight monazite grains from the ultramylonite AG17-J6b were analysed, for which concordant ²⁰⁶Pb/²³⁸U vs. ²⁰⁸Pb/²³²Th data (19) are between 596 Ma and 130 Ma (Fig. 7E). Three analyses at 500 Ma, 512 Ma and 596 Ma are obtained in a small-grained monazite from the matrix. Discordant data are located below and above the Concordia line and correspond to Th-poor (0.01 < Th/U < 0.88) and Th-rich monazites (1.99 < Th/U < 23.9), respectively. In the Th-rich monazites, ²⁰⁸Pb content mainly corresponds to radiogenic ²⁰⁸Pb, and common Pb incorporation generally results in an increase in the (common) ²⁰⁶Pb (analytical points above the Concordia). In the Th-poor monazites, radiogenic ²⁰⁸Pb is very low compared to ²⁰⁶Pb, and common Pb incorporation in these monazites results in the increase of the ²⁰⁸Pb (analytical points below the Concordia). It is thus possible to draw two Discordia lines for the two types of M2 monazites (Fig. 7E). Both Discordia independently (i.e., calculated with analyses obtained in distinct monazite grains) intercept the Concordia curve at the same age of 122 Ma.

The Tera–Wasserburg diagram clearly shows the increase in common Pb contamination in the M2 monazites (Fig. 7F). Concordant data spread between 168 Ma and 267 Ma with three inherited concordant data at 499, 508 Ma and 596 Ma. Four discordant data in the M2 monazites can be used to draw a Discordia line with a lower intercept age of 108 ± 17 Ma (MSWD = 1.8) and a higher intercept consistent with the common Pb isotopic composition at the same age (Supplementary Material S2).

5 Mica ⁴⁰Ar/³⁹Ar laser probe dating

Both in situ and step-heating ⁴⁰Ar/³⁹Ar experiments were conducted for this study on biotite and muscovite from seven gneiss samples showing variable degrees of mylonitic deformation (Tab. 1; Supplementary Materials 3 and 4). For step-heating analyses, a grain size between 150 µm and 400 µm was used and thin section observations show the dated micas are mainly clasts showing variable degrees of deformation-recrystallization preserved within the mylonitic foliation containing minute syn-kinematic muscovite and biotite. In situ ⁴⁰Ar/³⁹Ar analyses were carried out to compare the age of pre-kinematic and syn-kinematic micas.

Age spectra obtained on biotite and muscovite from these samples are shown in Figure 8. All the muscovites, regardless of the degree of mylonitisation, yield a narrow range of plateau ages from 103.4 ± 0.5 Ma to 106.7 ± 0.5 Ma for 80 to 100% of the ³⁹Ar released. Among these nine plateau ages, seven are very concordant at 104.7 ± 0.5 Ma. For these muscovites, despite a high radiogenic yield, the (³⁶Ar/⁴⁰Ar vs. ³⁹Ar/⁴⁰Ar) inverse isochron plots provide intercept ages that are consistent with the plateau ages, with no evidence of excess argon.

In five samples, which range from proto-mylonites to mylonites, biotites provide six age spectra. In the less deformed samples (samples AG17-J1, J2, J5a and b), the plateau ages scatter between 98.0 ± 0.8 Ma and 110.8 ± 0.7 Ma, with clear evidence of argon loss for two samples (AG17-J1 and AG17-J5a), related in the latter case to the presence of calcic impurities. The oldest biotite, with an age of 111 Ma (sample AG17-J5b), has a high, but imprecise, value of around 800 of the initial ⁴⁰Ar/³⁶Ar on the inverse isochron plot, suggesting that this mica has probably been contaminated by excess argon. Mylonite AG17-A2 provides two concordant plateau ages of 103.5 ± 0.5 Ma. ³⁸Ar_Cl is lacking in two samples (AG17-J5a and b) and present in the other three in a constant proportion with respect to the ³⁹Ar released during the step-heating procedure.

In situ ⁴⁰Ar/³⁹Ar ages obtained on 4 thin sections of proto- to ultra-mylonites (samples AG17-J2, J5b, J6a and A2) are summarized in Figure 8H, which shows cumulative probability plots of individual ages obtained for both large, partially recrystallized clastic and tiny syn-kinematic muscovite and biotite. The main observation is that texturally earlier micas have older ages than late micas in the foliation, as shown by the two age envelops drawn in Figure 8H. Large micas preserved within the foliation or in pressure shadows have ages mainly between 102 Ma and 112 Ma (n = 12) with a single age at 117 Ma, while the tiny micas in the foliation show a broader age range, falling into two groups at 94–102 Ma (n = 17) and 105–110 Ma (n = 7) respectively and a single age at 114 Ma. Neither step-heating nor in situ ⁴⁰Ar/³⁹Ar analyses provide any evidence of preserved Variscan ages.

Fig. 8 40Ar-39Ar laser step heating age spectra of muscovite and biotite. (A) Sample AG17-J1; (B) Sample AG17-J2; (C) Sample AG17-J5a; (D) Sample AG17-J5b; (E) Sample AG17-J6b; (F) Sample AG17-J7; (G) Sample AG17-A2; (H) In situ laser probe ages of pre- and syn-kinematic micas from variably deformed samples (AG17-J2, AG17-J5b, AG17-J6a and AG17-A2).

6 Interpretation of geochronological data

6.1 Significance of the different age groups

The oldest ages obtained from monazites between 500 ± 22 Ma and 596 ± 24 Ma, in samples AG17-J6b and AG17-J15b, undoubtedly represent a Cadomian inheritance as proposed by Fonteilles (1970) and Autran et al. (1970), and more recently by Tournaire Guille et al. (2019). In the M1 monazites included in garnet, two analytical points in samples AG17-A2 and AG17-J15B yield ²⁰⁶Pb/²³⁸U–²⁰⁸Pb/²³²Th concordant ages at 356 ± 9 Ma and 364 ± 9 Ma, which probably represent a mixing age between an early Variscan event in the mid-Ordovician and the main Variscan event during the late Carboniferous (Laumonier et al., 2008). Analyses carried out on these tiny monazite inclusions do not provide enough information for further interpretation.

In the literature, most of the ages obtained for the Agly Massif rocks and other North Pyrenean Massifs lie between 280 Ma and 320 Ma. They are commonly interpreted as reflecting both the influence of the Variscan HT-LP metamorphic event and the intrusion of magmatic rocks (Fonteilles, 1970; Respaut and Lancelot, 1983; de Saint Blanquat et al., 1990, de Saint Blanquat, 1993; Althoff et al., 1994; Olivier et al., 2004, 2008; Poujol et al., 2010; Odlum and Stockli, 2019, 2020; Tournaire Guille et al., 2019, Siron et al., 2020). Thus, we consider the ages between 280–320 Ma obtained in our samples as being related to the Variscan HT/LP metamorphic event. Interestingly, sample AG17-H1, which is a moderately deformed granite, yields only Variscan ages without any record of later M2 monazite growth (Figs. 6E and 6F).

The youngest concordant ages spread from around 280 Ma to 130 Ma, without showing any predominant population among the studied samples, despite the high number of analyses. Because no tectono-metamorphic event is known to have occurred during this period in this region, these ages can be interpreted as the result of i) lead loss by Variscan monazites inducing monazite rejuvenation, or ii) mixing between monazite domains of different ages during the LA-ICPMS analysis. Lead diffusion is extremely slow in monazite (Cherniak et al., 2004; Gardès et al., 2007) and volume diffusion is well known to be inefficient compared to dynamic dissolution-recrystallization processes in metamorphic and highly deformed rocks (Wawrzenitz et al., 2012; Didier et al., 2013, 2014). On the contrary, monazite rejuvenation in our samples is usually associated with common Pb contamination, which is clearly visible in the data points located above the Concordia curve in the Tera–Wasserburg diagrams. We thus interpret the ages between 280 Ma and 130 Ma as mixing ages. This mixing is due to the small size of the M2 monazite overgrowths or single grains (circa 5–20 μm) relative to the laser spot size (9 μm). The scattering reflects a variable amount of mixing between the M1 and M2 monazite domains. In this age range, the Concordia is very slightly curved and the Discordia line between the two mixing poles (i.e., Variscan ages and Cretaceous ages) merges with the Concordia curve. So, the analytical points appear to be concordant in the ²⁰⁶Pb/²³⁸U versus ²⁰⁸Pb/²³²Th diagrams.

The youngest ages obtained for our samples are Cretaceous (concordant ²⁰⁶Pb/²³⁸U–²⁰⁸Pb/²³²Th ages) and are recorded by the M2 monazites as well as by muscovites and biotites (⁴⁰Ar/³⁹Ar plateau ages). The most likely precursor phase for the Cretaceous M2 monazite is the Variscan M1 monazite since most of the M2 monazite grains crystallized as clusters surrounding M1 grains (Fig. 5D) or as overgrowth (Figs. 5C, 5E and 5F). In this context, dissolution-recrystallization processes in the presence of fluids and deformation can easily lead to the crystallization of new generations of monazites and/or micas with distinct compositions and ages. A regional-scale Cretaceous hydrothermal event has been well described and dated in the Agly Massif, as well as in the northern Pyrenees (Boulvais et al., 2007; Poujol et al., 2010; Fallourd et al., 2014; Boutin et al., 2016; Corre et al., 2018). Ages obtained for monazites and titanites from albitized samples are 110 ± 8 Ma and 98 ± 2 Ma, respectively, interpreted as the age of the albitization process in the Agly Massif (Boulvais et al., 2006; Poujol et al., 2010). Similarly, a Cretaceous ⁴⁰Ar/³⁹Ar age of 117.5 ± 0.4 Ma was obtained on a newly formed muscovite identified in an albitic matrix of the Hercynian Salvezines granite (Boulvais et al., 2007). These ages from non-mylonitized albitized samples are fully in line with those measured here in highly deformed non-albitized samples.

6.2 Deformation and fluids effects

Dissolution-recrystallization processes usually require fluids (Putnis, 2002; Putnis and John, 2010; Bosse and Villa, 2019). There is evidence of fluid-rock interaction in some of our samples. Garnet and locally feldspar display resorbed textures along rims or cracks, replaced by secondary biotite-muscovite-chlorite aggregates associated with small M2 monazites (Figs. 5C and 5D). To explain the high Y and HREE contents in M2 monazites, an external contribution is required in view of depleted Y content of hydrothermal monazite in albitized samples (Poujol et al., 2010). Resorption of garnet as observed in most of our samples is an efficient mechanism for the release and transport of Y and HREE in the presence of fluids. Likewise, Th-poor M2 monazites, as found in sample AG17-J6b, are well known in hydrothermal contexts (Grandhomme et al., 2016). All these observations point to evidence for synkinematic fluid-rock interaction in our samples. However, the chemical difference between hydrothermal monazites from albitized rocks and our monazites from mylonitic rocks suggests that the fluid-rock interaction documented in our samples is totally unconnected to the regional albitization event.

Thin section observations strengthen this assumption. Mylonite samples from the Agly Massif present microstructural evidence of conjugate top-to-the-NNE and top-to-the-SSW shear bands that developed in the presence of fluids. M2 monazite, as well as secondary muscovite and biotite (Fig. 5C), usually grow along these microstructures (Figs. 5 and 9). Sample AG17-J15b is characterized by ductile deformation at a relatively high temperature (around 500 °C) as well as by cataclastic deformation (around 200 °C) along the same deformation continuum (see the “Temperature conditions of the MPD” section). Figure 9 shows M1 monazite and K-feldspar, which are both involved in a shear band. As the temperature dropped, “dominoes” developed jointly in feldspar and monazite at about 450° C. These slip surfaces explain the apparent lengthening of the M1 monazite parallel to the foliation. BSE images show that the M1 monazite is partially recrystallized along these slip surfaces. Rims and recrystallized domains (i.e., M2 monazite) display concordant ²⁰⁶Pb/²³⁸U–²⁰⁸Pb/²³²Th ages of between 121 ± 4 Ma and 104 ± 3 Ma. The M1 monazite has recrystallized around its rim and along the shear planes (ramifications) resulting in its rejuvenation especially along the shear planes. The same feature is observed in sample AG17-A2: the M1 monazite has recrystallized in the garnet fractures, inducing the growth of a ramification with the youngest ²⁰⁸Pb/²³²Th ages at 116 ± 3 Ma, 112 ± 3 Ma and 111 ± 3 Ma (Fig. 5C). Finally, in sample AG17-J6b, M2 monazite clusters aligned along the main foliation yield ²⁰⁶Pb/²³⁸U ages of 124 ± 4 Ma and 127 ± 5 Ma (Fig. 5D). These results highlight the close bond between the Cretaceous ages recorded by the monazites and the structural features of the main mylonitic deformation. We can therefore conclude that M2 monazites crystallized due to deformation in the presence of fluids.

Step-heating and in situ ⁴⁰Ar/³⁹Ar mica ages in the range 94–107 Ma are in line with this interpretation. Dating shear zones in polymetamorphic basement rocks by the ⁴⁰Ar/³⁹Ar method is always a challenge due to several competing processes that affect the behaviour of argon in minerals during their isotopic resetting and closure. It has long been demonstrated that such behaviour is dependent not only on temperature but also on factors such as deformation, recrystallization, fluid-rock interaction, chemical and physical mineral properties, cooling rate and pressure (e.g., Dunlap et al., 1991; de Jong et al., 2001; Mulch et al., 2005; Di Vincenzo et al., 2006, 2016; Sanchez et al., 2011; Warren et al., 2012; Villa et al., 2014; Laurent et al., 2017; Abd Elmola et al., 2018; Oriolo et al., 2018; Bosse and Villa, 2019). As noted above, the stepwise heated micas in the Agly Massif are Variscan clasts preserved in the main foliation that experienced shearing and partial recrystallization, and the question remains as to whether the plateau ages of 98–105 Ma and in situ ages of 94–112 Ma reflect complete argon resetting or partial argon re-equilibration. Textural and mineral observations suggest that the temperature in the Agly basement at the beginning of deformation was as high as 500 °C and decreased progressively to below 300 °C by the end. Assuming that a high thermal gradient prevailed before shearing and therefore that argon was initially lost mainly by thermally activated volume diffusion in open system rocks, we used a diffusion loss model (Warren et al., 2012; Skipton et al., 2018) to evaluate the temperature required for full resetting of the muscovite and biotite clasts. For muscovite grains with radii between 200 µm and 50 μm, the temperature for total argon loss ranges from about 500 °C to 430 °C, with an uncertainty of about 50 °C due to the uncertainly on the activation energy of muscovite (63 ± 7 kcal/mol, Harrison et al., 2009). For biotite, this temperature decreases to 340–350 °C (± 25 °C) for a grain radius of 75–100 μm, which is consistent with its higher sensitivity to thermal events than muscovite. It is noteworthy that two in situ ages older than 112 Ma probably indicate that a small amount of inherited argon is still preserved in some mica grains despite a long thermal history coeval with deformation.

Therefore, these temperature estimates suggest that thermal conditions reached at the beginning of Cretaceous deformation in the Agly Massif were theoretically high enough to promote a total resetting of the ⁴⁰Ar/³⁹Ar chronometer in Variscan clasts of muscovite and biotite. This is consistent with the study of Odlum and Stockli (2019), which concludes that the temperature in the lower units of the Agly Massif was higher than 450 °C in the Aptian. Moreover, it is likely that this resetting was facilitated by continuous intra-grain deformation and recrystallization that resulted in the reduction of diffusion distances (e.g., Mulch and Cosca, 2004; Sanchez et al., 2011; Villa et al., 2014). The lack of significant age gradients in the age spectra (with the exception of two biotites that show evidence of recent argon loss) and the intra-sample and inter-sample reproducibility of plateau ages regardless of the deformation gradient tend to confirm that a high temperature gradient before and at the time of shearing played a major role in controlling the ⁴⁰Ar/³⁹Ar record of biotite and muscovite. Respective closure temperatures of these micas have been calculated using the diffusion parameters of Grove and Harrison (1996) and Harrison et al. (2009). For cooling rates in the range 10–20 °C/Ma, muscovite with a radius in the 50–200 μm range has a closure temperature between 380 °C and 430 °C while for biotite with a radius of 75–100 μm this temperature is close to 290–310 °C.

Based on temperature estimates from thin section analysis, this means that the 112–94 Ma mica ages can be interpreted as recording progressive recrystallization during continuous deformation from the moment that the ⁴⁰Ar/³⁹Ar chronometer was completely reset (step 1 of deformation) to closure temperatures for those micas. The later 94 Ma ages probably correspond to the end of syn-kinematic recrystallization of micas. These ages can be compared to those reported by Boulvais et al. (2007) for muscovites from strongly albitized Variscan granite and albitite from the Salvezines Massif, west of the Agly Massif. Magmatic muscovite displays a discordant saddle-shaped age spectrum ranging from 210 Ma to 240 Ma, suggesting that this mica experienced only partial resetting during albitization at about 400 °C. In contrast, the newly formed hydrothermal muscovite in the albitite has a well-defined plateau age of 117.5 ± 0.4 Ma, interpreted as dating its crystallization during albitization. This age matches the maximum mica ages reported in our study of mylonitic samples, which suggests that the beginning of deformation in the Agly middle/lower crust and Na-metasomatism in the upper crust were contemporary. Moreover, data reported by Boulvais et al. (2007) clearly show that primary muscovites from highly albitized rocks display only partially reset ⁴⁰Ar/³⁹Ar ages. In contrast, the mylonitic rocks dated in the Agly Massif show no evidence of metasomatism, which offers a strong argument that the Cretaceous mica ages cannot result only from a complete reset of the argon chronometer by interaction with fluids.

Fig. 9 (A) BSE image of monazite (sample AG17-J15b) with circles showing the location of the LA-ICPMS pits (9 µm) and their corresponding 208Pb/232Th ages at 2σ level (in red: concordant Cretaceous ages and in white: mixing ages). (B) A photomicrograph (top) and a diagram (bottom) showing the petro-structural environment of the analyzed monazite.

7 Variscan versus Cretaceous deformations

The geochronological results of this study contribute to the long-standing scientific debate concerning the Pyrenees as to whether the deformation and metamorphism in the Paleozoic units of the North Pyrenean Zone can be attributed entirely to the Variscan orogeny or whether they might be linked, at least partially, to Cretaceous extension (e.g., Passchier, 1984; de Saint Blanquat et al., 1986; Costa and Maluski, 1988; Paquet and Delay, 1989; Paquet and Mansy, 1991; de Saint Blanquat et al., 1990; Bouhallier et al., 1991; Monié et al., 1994; Olivier et al., 2008; Clerc et al., 2015; Odlum and Stockli, 2020; Vanardois et al., 2020). If part of the deformation in the LGU occurred during the Cretaceous rifting event, as suggested by geochronological data, can it be separated in the field from the earlier Variscan one?

A first observation is that Th-U/Pb Variscan ages are preserved in least deformed rocks, where they may even be the only ones recorded, as in sample AG17-H1 (Figs. 6E and 6F). Conversely, Cretaceous ages are systematically recorded in highly deformed zones associated with the strong N20E stretching lineation, indicating that the MPD in the LGU should be attributed to the rifting event.

In order to clarify the scientific debate, we compare our study with two recent, detailed ones that identified three successive deformations in the Agly Massif (Vanardois et al., 2020; Siron et al. 2020). The first one (D1) is a relictual fabric. This fabric, preserved in low-strain zones from the successive deformations, bears a mineral lineation oriented NW-SE. The second deformation (D2) is a flat-lying foliation that holds a conspicuous N10-20 trending stretching lineation. High strain zones from mylonites to ultramylonites can be observed in the whole structural sequence. Shear zones are parallel or slightly oblique to the main foliation and reveal a bulk top-to-the-north shearing with top-to-the-south conjugate shear zones. A third deformation (D3) is localized along a 4 km-wide E-W trending corridor, which crosscuts the whole massif. In this zone, the flat-lying D1/D2 foliation is progressively folded and rotated to the vertical. A vertical schistosity may develop parallel to the axial plane of the folds. Stretching lineations are very scarce and have variable plunge angles. Pinched Mesozoic slices like the Serre de Verges, Roquo-Roujo and the Serre de Cors are observed within the vertical gneiss series (Fig. 11 in Vanardois et al. 2020).

The authors consider that these three deformations (D1, D2 and D3) are of Variscan age. Neither the Cretaceous rifting nor the formation of the Pyrenees mountain chain had any imprint on the Agly Paleozoic Massif whose tectonic history concluded at the end of the Variscan orogeny (Vanardois et al., 2020).

Our field and thin section observations match this structural study closely. The relictual foliation of the D1 deformation has been observed in thin sections where it is folded and transposed into a new flat-lying foliation (Fig. 2E). There is no doubt that this early deformation, which is almost erased, is Variscan. The D2 deformation which displays flat-lying foliation, a strong NNE-SSW trending stretching lineation associated with high strain zones of opposite top-to-the-NNE and top-to-the-SSW shearing corresponds to the MPD that we observed in the field, analyzed in thin sections and dated as Cretaceous by two geochronological methods. The discrepancy with the study of Vanardois et al. (2020) is that we did not identify a bulk top-to-the-NNE shearing, but instead found equal top-to-the-NNE and top-to-the-SSW shearing. For instance, a 100 m-thick top-to-the-SSW shear zone in a highly deformed orthogneiss can be observed south of Planèzes in the Canorgues area, which shows an impressive strain gradient from bottom (AG17-A2, ultramylonites) to top (AG17-A3, moderately deformed).

In a companion paper, the same group of researchers carried out a study in which they determined the timing of the metamorphism in the Agly Massif using in situ monazite U-Th-Pb dating (Siron et al., 2020). The main findings concerned the granulite metamorphism dated at ca. 305–306 Ma while the D2 deformation of Vanardois et al. (2020) was dated at ca. 296–300 Ma based on ages in monazite overgrowths from a highly mylonitized kinzigite (sample 11Ag03 in Siron et al., 2020). The authors' conclusions focus on Variscan ages measured in pristine samples with little or no deformation. However, the authors note that “in all the analyzed samples”, “a spread in age along the concordia” is observed giving ages up to 146 Ma (Tab. 4 and Fig. 11 in Siron et al. 2020). The more deformed the sample, the younger the ages. Finally, they “speculate” that these ages are “the result of localized fluid percolation during Cretaceous event (ca. 100 Ma)”, despite the absence of overprint in their samples. This speculation is not commented on further in their discussion.

These results are in perfect agreement with the results presented here, considering that the samples studied by Siron et al. (2020) are not or only slightly deformed. The age dispersion observed in their samples is similar to that obtained in the least deformed sample AG17-H1 (Figs. 6E and 6F). The so-called “highly mylonitized sample” (11Ag03 in Fig. 9 in Siron et al., 2020) is in fact a very heterogeneously deformed sample. Narrow high-strain shear bands criss-cross the thin section but the foliation is moderately deformed while garnets are completely undeformed (11AG03 in Fig. 9 in Siron et al., 2020). The five monazites analyzed all come from outside the narrow highly strained zones. Accordingly, these monazites yield the age of the Variscan foliation as our sample AG17-H1, and not the age of the Cretaceous shearing deformation that we obtained in the highly sheared mylonites where garnets are strongly deformed (Fig. 10A). Again, this illustrates that the Cretaceous deformation reactivated the former Variscan foliation.

The age of the main flat-lying surface (S2 for Vanardois et al. 2020, MPD this study) is not unequivocal. This flat-lying surface bears mineralogical assemblages formed during the HT/LP Variscan event. From a metamorphic point of view, this flat-lying surface is Variscan, even if HT/LP assemblages were partially replaced later by a syn-kinematic assemblage including biotite and muscovite.

However, when highly sheared, geochronological data indicates that strain, including stretching lineations and kinematic indicators, is Cretaceous. From a kinematical point of view, this surface is Cretaceous. We would like to stress that metamorphism and strain should not be confused. The Variscan foliation has been reused and reworked during the rifting event (and sometimes isoclinaly folded and transposed into a new flat-lying surface, Fig. 2E). Reusing the Variscan foliation, there is no reason why the colder Cretaceous deformation should remove the Variscan metamorphic assemblages.

In this interpretation, a question remains open. According to this geochronological study, the N20E stretching lineation in the Agly Massif is associated with a Cretaceous shearing deformation. This is in good agreement with the N20E stretching lineation measured in the surrounding Mesozoic series (Bas Agly and Boucheville synclines). It also fits the opening direction of the Cretaceous Rift. However, according to Vanardois et al. (2020), preferential alignments of cordierites and andalusites define a N20E direction. Further studies are needed to determine whether these minerals were rotated during the Cretaceous event or whether they reveal the direction of the Variscan stretching lineation, which in that case would have the same orientation as the Cretaceous.

As we have avoided working in the steeply dipping E-W corridor, we have no relevant information on the D3 deformation. We have observed it east of Rasiguères in the Devonian calcschists where a first schistosity is intensively deformed into folds forming a new cleavage by transposition of the first schistosity into vertical axial planes (Fig. 10B). No evidence can be found to determine if the first schistosity in this area is an S1, S2 or an S1-S2 surface. However, pinched Mesozoic slices within the E-W vertical corridor (D3 deformation) strongly suggest that it formed during the Pyrenean orogeny. Also, its orientation similar to that of the Agly and Boucheville synclines argues in favor of a Pyrenean event.

To conclude, we observed the same structural relationship between three successive deformations in the field as Vanardois et al. (2020) and our geochronological study strongly suggests that the D2 deformation (i.e., the MPD of this study) is related to the Cretaceous event. Accordingly, the D3 compressional deformation with pinched Mesozoic slices located in an E-W corridor where previous schistosities were rotated into a vertical position should be ascribed to the formation of the Pyrenean chain.

Fig. 10 (A) Strongly deformed garnet in an ultramylonite of the MPD deformation in AG17-A2. See also fig. 2D and compare with undeformed garnets of Figure 9 in Siron et al. (2020). Picture 6 (B) refers to Figure 1 where the location of the photos is shown (orange circle): post MPD deformation: folded Devonian calcschist (east of Rasiguères).

Fig. 11 Interpretative PTt path for the LGU unit. Note the contrasting geotherms and the diachronism of the temperature peak between the LGU and the sediments of the two sedimentary basins (Explanations in the text).

8 PTt path of the LGU

The crust in the Agly Massif is considered to be 20 km thick at the end of the Variscan orogeny with PT conditions in the basement of 600/750 °C at a depth of about 10/15 km (Guittard and Raquin, 1958; Fonteilles, 1970; Vielzeuf and Kornprobst, 1984; Delay, 1989; Delay and Paquet, 1989; Paquet and Mansy, 1991; Olivier et al., 2004; Cochelin et al., 2017; Siron et al., 2020). Starting from this HT-LP Variscan metamorphism at 300 Ma (point A in Fig. 11), we attempt to reconstruct the PTt path of the LGU for a point buried at a depth of 15 km at that time.

It is worth noting that there was a gap of nearly 200 Ma between the HT-LP Variscan metamorphism (300 Ma) and the very start of the Cretaceous rifting (around 120 Ma). Cooling of the lithosphere following extensional regimes lasts a few tens of millions years, and in all cases is completed within 100 Ma (e.g., McKenzie, 1978; Jarvis and McKenzie, 1980; Allen and Allen, 2018). Complete thermal re-equilibration over 200 Ma is likely for the Pyrenean lithosphere, which had a thin crust (20 km) at the end of the Variscan orogeny. No surface erosion is needed to complete thermal re-equilibration, just a shift of the thermal boundary at the base of the lithosphere to a deeper level. For this reason, we suggest that the crust returned to an undisturbed geotherm typical of anorogenic continental crust long before the start of the Cretaceous rifting (point A’ in Fig. 11). Accordingly, cooling without erosion causes point A to move to the undisturbed geotherm (point A’ in Fig. 11).

However, according to Odlum and Stockli (2019), the temperature of the middle and lower crust in the Agly Massif was already higher than 450 °C in the Aptian at the beginning of the rifting. This is in agreement with our temperature estimate (about 500 °C) deduced from thin section analysis for the beginning of the Cretaceous deformation. Exhumation and cooling of the middle-lower crust (i.e., the LGU) started in the early Aptian (Odlum and Stockli, 2019). This is consistent with the oldest age of this study (127 Ma) which points to a Barremian age for the very start of the rifting event. Exhumation and cooling continued through the Albian and Cenomanian (Odlum and Stockli, 2019). The Agly Massif then experienced rapid cooling to 200–250 °C by 90 Ma (Ternois et al., 2019). Again, this is consistent with the youngest ages of this study (94 Ma), which probably represent the end of syn-kinematic recrystallisation of micas in the LGU. This gives a range for the duration of the extension in the Agly Massif from about 130 Ma to 90 Ma.

The Cretaceous rifting would have been responsible for the thinning of the Agly Massif from 20 km to 8 km thick (Clerc and Lagabrielle, 2014; Clerc et al., 2015; Ternois et al., 2019) but the variation in the thinning rate from 130 Ma to 90 Ma is unknown. For this reason, we represent the PTt evolution in the LGU during the Cretaceous rifting by an approximate anti-clockwise path showing firstly reheating of unknown origin to 500 °C in the Barremian, then thinning, cooling and exhumation from 130 Ma to about 90 Ma (Fig. 11). While this anti-clockwise path is relatively well constrained by age-temperature relationships, the age-depth pattern remains less well defined.

This PTt evolution in the LGU is in marked contrast with the geotherm deduced from the maximum temperature of the HT-LP Cretaceous metamorphism recorded in the Boucheville and Agly basins (Fig. 11). Also, the thermal peak is reached later (95 Ma) (e.g., Clerc et al., 2015) compared to the Agly Massif, which is a strong constraint in the proposition of a crustal-scale tectonic model.

9 Crustal-scale evolution

This study suggests that the Cretaceous thinning in the Agly Massif occurred through numerous more or less parallel but opposite shear zones across the entire sequence of the LGU. This does not preclude a significant part of Cretaceous thinning also occurring along the Latour-Caladroy detachment fault observed at the boundary between the LGU and the USU. It is known that mechanical decoupling between the upper and lower crust in the Agly Massif took place by activation of this fault, along which 5.5 km may have been subtracted (e.g., Delay, 1990; Bouhallier et al., 1991; Siron et al., 2020).

The Agly Massif is limited to the north and south by two large shear zones at the contact with the Boucheville and Agly Mesozoic series. To the north, this is shown by top-to-the-NNE shear features in mylonites observed in the lower part of the Mesozoic sediments next to the Agly Massif (Vauchez et al., 2013). To the south, the highly deformed Mesozoic sediments along the contact have been interpreted as being the result of a top-to-the-SSW extensional fault zone, now tilted vertically, or even inversed at some places, by the later compressional Pyrenean deformation (Chelalou et al., 2016).

A N20E oriented stretching lineation is observed in the Mesozoic units of the Agly and Boucheville basins (Vauchez et al., 2013; Chelalou et al., 2016). According to these authors, this stretching lineation is contemporaneous with the Cretaceous thermal event observed in the Pyrenean Mesozoic units (110–90 Ma). Combined with our results, this indicates that both the basement and the Mesozoic cover underwent the same extensional deformation in the Cretaceous.

It follows that the Agly Massif appears to be a ductile boudin bounded by two large-scale detachment faults, which led the basement to complete rupture below the Boucheville basin to the south and thin to a narrow neck below the Agly basin to the north, as documented by the pressure/temperature conditions recorded in the sedimentary units (see references therein).

The two cross-sections in Figure 12 enhance the diachronism of the Cretaceous metamorphic evolution between the Agly Massif and the two bounding sedimentary basins. From the start of the extension (i.e., in the late Barremian or early Aptian), there was active cooling in the Agly Massif (Odlum and Stockli, 2019 and this study) while the formation of necking zones on both sides of the massif increased the temperature in sedimentary units (Fig. 12A). In the Cenomanian, rupture of the basement below the Boucheville basin and its severe thinning below the Agly basin were contemporaneous with peak metamorphism in the sedimentary units at about 95 Ma, while the Agly Massif, which lay between the two, cooled rapidly to below the ductile/brittle transition (Ternois et al., 2019) (Fig. 12B).

Fig. 12 Idealized NS cross-sections of the basement crust deformation in the Agly area during the Cretaceous rifting. (A) Aptian. The basement crust deformed as a pinch and swell structure, leading to an isolated boudin of the Agly Massif bounded by two detachment faults. The Agly Massif was cooling from the start of the rifting. (B) Albo-Cenomanian. In the Cenomanian (95 Ma), the basement crust reached the point of rupture below the Boucheville basin, which caused a late HT-LP metamorphism in the Mesozoic sediments (up to 580 °C at a depth of 7 km). In the Agly basin, the thinning of the basement crust was less extreme, but enough to be responsible for a late temperature overprint (95 Ma) of up to 390 °C in the Mesozoic sediments. NPF: transform fault separating the Iberian microplate from Europe (= future North Pyrenean Fault); LCDF: Latour-Caladroy Detachment Fault (explanations in the text).

10 Inherited heat flow prior to the start of the rifting

The modalities of the extension point to a lithospheric stratification with strong viscosity contrasts generating large scale pinch and swell structures in the basement of the North Pyrenean Zone. For this process to occur, the basement has to be sandwiched between weaker strata. As the extension progresses, these low strength horizons may flow in a ductile manner to fill the gaps created by necking zones in the basement crust (Fig. 12).

The thick layer of very weak Triassic sediments (up to 2700 m), just above the basement, is an excellent candidate, even at low temperatures, for creating a complete decoupling between the sedimentary upper crust and the underlying basement crust, as already proposed by some authors (e.g., Lagabrielle et al., 2010, 2020; Asti et al., 2019; Duretz et al., 2019).

It is less clear what the situation at the base of the crust is. Where the mantle is exposed in the North Pyrenean Zone, the Variscan granulitic lower crust is generally not present or restricted to very small, highly disrupted outcrops (e.g., Corre et al., 2016; Asti et al. 2019). This indicates its (nearly) total disappearance due to extreme thinning. Furthermore, the final geometry in the necking zones, showing that the mantle and the Trias (and other pre-rift sediments) are in contact, implies that the mantle also flowed and behaved as a mechanically weak horizon. The fact that the mylonitic deformation in the mantle started at high temperatures (950 °C in the western Pyrenees and 800 °C in the Turon de la Técouère Massif in the eastern Pyrenees) and low pressures (5–10 kbars) (Fabries et al., 1991; Vissers et al., 1997; Newman et al., 1999) strengthens this point of view. Likewise, the lack of melting or refertilization of the mantle exhumed during the Cretaceous rifting also confirms that it comes from shallow depths (Picazo et al., 2016). These high temperatures at shallow levels (20–40 km) reveal that the uppermost lithospheric mantle was already ductile at the start of the rifting episode. It follows that the viscosity contrast between the Variscan granulitic lower crust and the uppermost mantle must have been low or even nonexistent during the rifting.

Such stratification of the lithosphere is not consistent with an undisturbed geothermal gradient at the beginning of deformation. Considering that the continental crust in the Agly Massif was 20 km thick at the end of the Variscan orogeny, and therefore at the onset of the Cretaceous rifting, an average geothermal gradient in a thin crust would be associated with a strength profile involving a brittle upper crust and a brittle uppermost mantle. In between, the middle crust would have brittle/ductile transition conditions. For a thin crust (20 km), only its very base would become ductile in an anorogenic environment. This does not reflect the way the crust and the uppermost mantle deformed from the very start of the rifting. The middle crust underwent ductile boudinage while the thick Triassic evaporites and the uppermost mantle flowed in a very ductile manner to come into contact with each other. Such a strength profile for the lithosphere is only compatible with a high geothermal gradient

Indeed, the process of crustal lateral extraction requires a high geothermal gradient (See review in Clerc et al., 2018). Elevated heat flow within the lithosphere is an imperative condition for promoting viscosity contrasts high enough to trigger crustal boudinage of the middle crust. Consistent with the relatively high temperatures in the late Barremian, > 450 °C or > 500 °C (Odlum and Stockli, 2019 and this study, respectively), high heat flow may result from a blanketing effect of the syn-rift sediments above the basement or from an inherited high geothermal gradient (Clerc et al., 2018).

A blanketing effect (i.e., thermal insulation of the basement below thick sedimentary series) can be precluded as the extension is already well under way in the Aptian, which is associated with thin sedimentation (< 500 m, Berger et al. 1993) that is insufficient to produce such an effect. Only the very thick Albo-Cenomanian sedimentation (up to 5000 m) could have produced a blanketing effect, but at that time the Agly Massif was cooling rapidly. An inherited heat flow is an alternative hypothesis. This may have been due to infiltration of hot fluids and/or upward heat transfer along the vertical lithospheric limit between the Iberian and European plates. However, the investigation of this hypothesis lies beyond the scope of this study.

11 Conclusion

This work is a new piece in the puzzle that, progressively, reveals the major tectonic significance of the Pyrenean continental rifting during the Cretaceous. The main results of this study of the Low Gneiss unit (LGU) in the Agly Massif can be summarized as follows:

1. The Variscan foliation bearing the HT/LP metamorphic assemblages was reactivated and reworked in the Cretaceous by a highly shearing deformation, which, together with the formation of new synkinematic minerals, was associated with an omnipresent N20E stretching lineation and numerous kinematic indicators.

2. This Cretaceous deformation developed from about 127 Ma to 94 Ma as indicated by Th-U/Pb of monazite and ⁴⁰Ar/³⁹Ar dating of mica retrieved from proto-to ultramylonites.

3. Cretaceous extension took place by inhomogeneous strain within the entire LGU unit. Coeval shear zones top-to-the-NNE and top-to-the-SSW, which flattened with time to become more or less parallel, resulted in lengthening of the whole massif in a NNE-SSW direction.

4. The Agly Massif is a ductile boudin limited by two necking zones, the Boucheville basin to the south and the Agly basin to the north. The diachronism in the metamorphic evolution between the Agly Massif and the sediments of both basins is explained by cooling of the Agly Massif from the beginning of the rifting while synchronous necking of the basement below the basins increased the temperature in those sediments due to mantle exhumation.

5. Ductile boudinage of the basement in an already thin crust implies an inherited heat flow in the Agly Massif prior to the start of extension. This is not due to a blanketing effect but rather to upward heat transfer along the lithospheric fault separating the Iberia and European plates.

Acknowledgements

This paper is a Laboratory of Excellence ClerVolc contribution with the number 467. We acknowledge Yves Lagabrielle for the fruitful discussions on the Cretaceous Pyrenean Rift. Authors would like to thank the editor Laurent Jolivet, the associated editor Romain Augier and two anonymous reviewers for their helpful reviews and remarks that improved the quality of our manuscript.

References

Supplementary Material

Supplementary Material S1. Analytical methods.

Supplementary Material S2. U-Th-Pb dating.

Supplementary Material S3. Argon in situ.

Supplementary Material S4. Argon step heating.

All Tables

All Figures

Fig. 1 (A) Simplified geological map of the North Pyrenean Zone (NPZ) and location of the Agly Massif in the Pyrenean belt. (B) Geological map (modified from Fonteilles et al., 1993 and Berger et al., 1993) showing the location of the Agly Variscan Massif between the Saint Paul de Fenouillet and Bas Agly sedimentary basins to the north and the Boucheville sedimentary basin to the south. The North Pyrenean Fault (NPF) separates the Northern Pyrenean Zone from the Axial Zone. Stretching lineations from Delay (1989), Bouhallier et al. (1991), Choukroune et al., (1990), Vauchez et al., (2013) and this study.

In the text

	Fig. 2 Pictures A, B, C and D refers to the Figure 1 where the location of the photos is shown (orange circle). (A) MPD flat-lying mylonitic shear zone in orthogneiss (east of Mas Chiffre-Rnes). (B, C and D) Kinematic indicators of the MPD (St-Martin/Cuxous area). (E) Thin-sections showing the folding of the relictual fabric (RF) transposed onto the MPD foliation.
In the text

Fig. 3 Thin sections showing (A) dominant top-to-the-SSW senses of shear (SC structures) at the upper half of the thin section and dominant top-to-the-NNE senses of shear at the lower half interpreted as the SSW ending of a microscopic lens, (B and C) symmetric MPD lenses around resistant feldspars, (D) detail of a late semi-brittle normal fracture crosscutting the whole thin section, (E) MPD domino-type bookshelf feldspars indicating parallel but opposite senses of shear. The thin section consists of parallel bands indicating alternating top-NNE and top-SSW shear senses. In each parallel band, bookshelf feldspars are consistent with the flow revealed by quartz fabrics, pressure shadows and shear bands, (F) outcrop showing centimetre-scale MPD lenses (Saint Martin area). Picture E refers to the Figure 1 where the location of the photos is shown (orange circle).

In the text

Fig. 4 Photomicrographs of quartzofeldspathic gneisses with biotites and garnets. (A) Almond-shaped K-feldspar surrounded by fine-grained aggregates of biotites and recrystallized quartz ribbons. (B) Almond-shaped K-feldspar later deformed into a domino-type. (C) Domino-type feldspar later partially ground at its tails. (D) Almond-shaped K-feldspar later deformed into a domino-type and finally completely disrupted in a cataclastic manner at its tails. The feldspar is surrounded by fine-grained quartz ribbons. (E) K-feldspar crosscut by conjugate Coulomb-type fractures dipping at 60 degrees. (F) Fractured garnet crosscut by conjugate Coulomb-type fractures dipping at 60 degrees.

In the text

Fig. 5 BSE images of monazite (samples AG17-A2, AG17-J15b and AG17-J6b). Monazites are found as inclusions within large garnets (A) or found in the matrix (B, C and D) in contact with biotites (D), K-feldspars (B) or garnets (C). Large monazite grains are surrounded by small grains or recrystallized branches (C) or can be clearly deformed (B). Small monazite grains can also form clusters aligned along the foliation plane (D). La (left) and Y (right) X-ray maps of monazites (E) in matrix and (F) at the edge of garnet in sample AG17-J15a. Circles show the location of the LA-ICPMS pits (9 μm) and their corresponding 208Pb/232Th ages. Circles show the location of the LA-ICPMS pits (9 μm) and their corresponding 208Pb/232Th ages (In red: concordant Cretaceous age; in blue: mixing age; in black: Variscan age).

In the text

Fig. 6 (A) 208Pb/232Th versus 206Pb/238U diagram for monazites in sample AG17-A2. Dashed ellipses correspond to discordant data. Dark ellipses correspond to concordant ages affected (or not) by mixing between M1 and M2 monazite domains. White ellipses correspond to M2 analyses affected by common Pb contamination and mixing with M1 domain. Grey ellipses correspond to analyses used for the Discordia calculation. (B) Tera–Wasserburg diagram for the same sample. (C) 208Pb/232Th versus 206Pb/238U diagram for monazite in sample AG17-A3. (D) Tera–Wasserburg diagram for the same sample. (E) 208Pb/232Th versus 206Pb/238U diagram for monazite in sample AG17-H1. (F) Concordia diagram for the same sample.

In the text

Fig. 7 (A) 208Pb/232Th vs. 206Pb/238U diagram for monazites in sample AG17-J15a. Dashed ellipses correspond to discordant data. Dark ellipses correspond to concordant ages affected (or not) by mixing between M1 and M2 monazite domains. White ellipses correspond to M2 analyses affected by common Pb contamination and mixing with the M1 domain. Grey ellipses correspond to analyses used for the Discordia calculation. (B) Tera–Wasserburg diagram for the same sample. (C) 208Pb/232Th vs. 206Pb/238U diagram for monazites (sample AG17-J15b) and (D) associated Tera–Wasserburg diagram. (E) 208Pb/232Th versus 206Pb/238U diagram for monazites in sample AG17-J6b. Dashed ellipses correspond to discordant data. Grey ellipses correspond to analyses used for the Discordia calculation. This sample is characterized by two populations of chemically distinct monazites, which allow two Discordia lines to be drawn. See the text for detailed explanations. (F) Tera–Wasserburg diagram for the same sample.

In the text

	Fig. 8 40Ar-39Ar laser step heating age spectra of muscovite and biotite. (A) Sample AG17-J1; (B) Sample AG17-J2; (C) Sample AG17-J5a; (D) Sample AG17-J5b; (E) Sample AG17-J6b; (F) Sample AG17-J7; (G) Sample AG17-A2; (H) In situ laser probe ages of pre- and syn-kinematic micas from variably deformed samples (AG17-J2, AG17-J5b, AG17-J6a and AG17-A2).
In the text

	Fig. 9 (A) BSE image of monazite (sample AG17-J15b) with circles showing the location of the LA-ICPMS pits (9 µm) and their corresponding 208Pb/232Th ages at 2σ level (in red: concordant Cretaceous ages and in white: mixing ages). (B) A photomicrograph (top) and a diagram (bottom) showing the petro-structural environment of the analyzed monazite.
In the text

	Fig. 10 (A) Strongly deformed garnet in an ultramylonite of the MPD deformation in AG17-A2. See also fig. 2D and compare with undeformed garnets of Figure 9 in Siron et al. (2020). Picture 6 (B) refers to Figure 1 where the location of the photos is shown (orange circle): post MPD deformation: folded Devonian calcschist (east of Rasiguères).
In the text

	Fig. 11 Interpretative PTt path for the LGU unit. Note the contrasting geotherms and the diachronism of the temperature peak between the LGU and the sediments of the two sedimentary basins (Explanations in the text).
In the text

Fig. 12 Idealized NS cross-sections of the basement crust deformation in the Agly area during the Cretaceous rifting. (A) Aptian. The basement crust deformed as a pinch and swell structure, leading to an isolated boudin of the Agly Massif bounded by two detachment faults. The Agly Massif was cooling from the start of the rifting. (B) Albo-Cenomanian. In the Cenomanian (95 Ma), the basement crust reached the point of rupture below the Boucheville basin, which caused a late HT-LP metamorphism in the Mesozoic sediments (up to 580 °C at a depth of 7 km). In the Agly basin, the thinning of the basement crust was less extreme, but enough to be responsible for a late temperature overprint (95 Ma) of up to 390 °C in the Mesozoic sediments. NPF: transform fault separating the Iberian microplate from Europe (= future North Pyrenean Fault); LCDF: Latour-Caladroy Detachment Fault (explanations in the text).

In the text