Pre‐Subduction Architecture Controls Coherent Underplating During Subduction and Exhumation (Nevado‐Filábride Complex, Southern Spain)

The interplay between structural and metamorphic processes operating along the deep plate interface in subduction zones remains elusive as much of the geologic record is recycled into the mantle. In some cases, metamorphosed subducted rocks are underplated and exhumed to the surface, providing critical constraints on structural processes and the rheological evolution of subduction interfaces at convergent margins. One such exhumed high‐pressure/low‐temperature subduction complex is the Cenozoic Nevado‐Filábride Complex (NFC) in Southern Spain. This study presents new data from the NFC that elucidate the syn‐metamorphic deformation, stacking, and underplating of continental slivers along the subduction interface. The structurally lowest NFC dominantly comprises lithologically monotonous Paleozoic metamorphic basement rocks recorded by apatite U‐Pb ages and shows no evidence for large‐scale internal duplications suggesting it behaved as a coherent basement succession during subduction. In contrast, structurally higher levels of the NFC are characterized by the stacking of older‐on younger coherent slices and distinctly different metamorphic ages. These relationships document syn‐subduction structural repetitions and tectonic stacking of imbricate thin slivers (∼100s m) during subduction underplating. Structurally higher levels of the NFC exhibit both Eocene and Miocene metamorphic zircon rims and apatite ages, along with microstructures indicative of relatively higher temperature metamorphism. Large‐scale underplating and antiformal stacking of slivers in the subduction channel can provide buoyancy forces to underplate and assist exhumation. We demonstrate that the presubduction stratigraphic architecture is a key control on the style and timing of deformation and metamorphism, facilitating coherent subduction underplating.

Reconstructing the presubduction crustal anatomy and stratigraphy from the initial deposition and characterizing how they are subsequently modified is critical for understanding the nature, style, and spatial scales of deformation during subduction and the coupling between the downgoing and overriding plates. Rocks that end up being subducted and underplated likely experience multiple stress regimes and deformation styles due to changing tectonic conditions during subduction and exhumation. Hence, most of the presubduction geologic record is typically overprinted and initial stratigraphic characteristics have been annealed. However, as our knowledge of these regimes improves, there have been recent discoveries and methodological advances that allow us to reconstruct the presubduction configurations in the exhumed rock record and provide insights into subduction processes. A few examples include preserved rift-related structures in the Alps and Pyrenees (Jammes et al., 2009) as well as in the Betic Cordillera (Martin-Rojas et al., 2009). Studies from the subduction complexes in the Cyclades, Greece (Kotowski et al., 2022;Poulaki et al., 2019;Seman, 2016) show that the presubduction architecture of these exhumed terranes can be discerned. Discriminating presubduction versus syn-subduction processes is essential to better understand the spatial and temporal scales of fault and shear zone motions, rheological conditions, and the specific pressure and temperature conditions during subduction zone metamorphism. Diffraction (EBSD) analyses and field observations, our data show that the initial stratigraphy of the NFC was rearranged during subduction, underplating, and exhumation throughout the Cenozoic and provide precise timing constraints on the record of metamorphism and initiation of the Betic subduction in the Western Mediterranean.

Mediterranean Subduction Complexes
Subduction zones in the Mediterranean region have a unique architecture that differentiates them from other subduction zones across the world. The major tectonic events of the Mediterranean include convergence and mountain building during the Carboniferous Variscan orogeny followed by rifting during the Triassic/Jurassic and subsequent opening of domains of the Neo-Tethys Ocean (Asti et al., 2019;Pedrera et al., 2020). In the Late Cretaceous, convergence between Africa and Europe, due to the anticlockwise rotation of Africa, triggered the Alpine orogeny with segmented subduction systems initiating from east to west (e.g., Dewey et al., 1989;Faccenna et al., 2004;Jolivet et al., 2003). Isolated oceanic and rifted continental blocks were subducted, underplated, and later exhumed during trench retreat (e.g., Brun & Faccenna, 2008;Lonergan and White, 1997;Scott & Lister, 1992). Due to deformation associated with collision and complicated outcrop geometries of the subduction complex rocks in alpine settings, as well as the disappearance of numerous micro-oceanic domains, the rapid evolution of the region makes it challenging to use traditional methods to reconstruct the geologic history of the Mediterranean.
One of the most well-exposed subduction complexes in the western Mediterranean is the NFC. The NFC contains the geologic record of the formation and breakup of Pangea and underwent subsequent subduction and exhumation in the Cenozoic. At present, the NFC outcrops along major folded elongated extensional domes in the Betic Cordillera of southern Spain . The Betic Cordillera forms the northern component of the arcuate Betic-Rif orogenic belt, whereas the southern counterpart comprises the Rif Mountain belt in northern Morocco. The general architecture of the Betic-Rif system was formed with the subduction of oceanic and continental crust due to convergence between Africa and Eurasia, followed by westward migration of the trench and coeval collision of the Alboran domain with the African and South Iberian paleomargins (e.g., Balanyá et al., 1997;Booth-Rea et al., 2005. Between the Betic-Rif Cordillera lies the Alboran upper plate domain in the westernmost Mediterranean ocean, which largely comprises of the Alpujarride Complex and has also experienced metamorphism and crustal shortening until the Early Miocene (e.g., Esteban et al., 2011;Platt et al., 1998). Offshore regions of the Alboran basin in the Betics and Rif were below sea level since the Early to Late Miocene (e.g., de la Peña et al., 2021;Rodríguez-Fernández et al., 2011).
Subduction-related volcanism was sparse throughout Cenozoic convergence but was active from ∼18 to 6 Ma and is currently exposed in southeastern Spain, providing evidence of the subduction geometry in the Miocene (Duggen et al., 2003(Duggen et al., , 2004(Duggen et al., , 2008Varas-Reus et al., 2017). The tectonic evolution of this system is complex because of the exposure of subduction-related terranes, both subcontinental and oceanic mantle sequences, and a strong component of oblique and transcurrent tectonics. Additionally, there is evidence for a significant amount of trench retreat that drove subsequent extension in the overriding plate (Faccenna et al., 2004;Gutscher, 2012;Lonergan & White, 1997) and lateral slab tearing at the edges of the system along the Betics and Rif (e.g., Capella et al., 2020;de Lis Mancilla et al., 2015;García-Castellanos & Villaseñor, 2011;Levander et al., 2014), leading to the present-day near-vertical position of the slab beneath the Alboran Sea (Bezada et al., 2013).

The Structural and Stratigraphic Rearrangement of the Nevado-Filábride Complex
The geologic terranes of the Betic Cordillera are grouped into Internal and External zones. The External zones contain Mesozoic to Miocene sedimentary rocks deposited along the Iberian continental paleomargin. While the Internal Betics have been divided into three main groups from top to bottom: the unmetamorphosed Maláguide, and the metamorphosed Alpujárride and NFC (Balanyá & García-Dueñas, 1987). Initially, they were all considered to be part of the Alboran domain (i.e., the allochthonous terrane translated between Iberia and Africa; Balanyá & García-Dueñas, 1987); however, more recent studies have suggested the NFC formed along the southern margin of the Iberian Peninsula (e.g., Booth-Rea et al., 2015;Gómez-Pugnaire et al., 2004, 2012Jabaloy-Sánchez et al., 2018, 2021Kirchner et al., 2016;López Sánchez-Vizcaíno et al., 2001;Platt et al., 2006;Poulaki & Stockli, 2022). are at least three different tectonic units separated by metamorphic grade with shear zones marking their contacts: the basal Ragua unit (Gómez-Pugnaire & Franz, 1988), Calar Alto, which includes the Tahal formation, and the upper Bédar-Macael unit. Another well-accepted subdivision considers two main units: a lower homogeneous Veleta unit and an upper, more heterogeneous Mulhacén succession, which includes metamorphosed mafic and ultramafic rocks (Puga et al., 2000(Puga et al., , 2002. Detrital zircon studies have shown that the Veleta unit is Carboniferous and older in age , 2021Poulaki & Stockli, 2022;Santamaría-López & Sanz de Galdeano, 2018). Carboniferous ages of the Veleta unit are also supported by conodont fossils (Rodríguez-Cañero et al., 2018). Additionally, Laborda-López et al. (2015) identified Early Devonian fossil assemblages in graphite marbles intercalated in the graphite schist of the Veleta unit. In contrast, detrital zircon data show that the overlying Mulhacén succession is Late Carboniferous/Early Permian to Early Jurassic , 2021Poulaki & Stockli, 2022). More recently, Poulaki & Stockli (2022) performed abundant detrital zircon dating and found that the NFC is Paleozoic to Early Jurassic metasedimentary sequence preserving the record of sedimentation from the Variscan orogeny to Neo-Tethys opening. The largest part of the Tahal formation dates to the Permian-Triassic/Early Jurassic (Poulaki & Stockli, 2022) with Jurassic basic intrusions (Puga et al., 2005). The Bédar-Macael subunit represents the structurally highest parts of the NFC. The Bédar-Macael includes various types of schist primarily Permian in age (Poulaki & Stockli, 2022), thick carbonate units that are generally assumed to be Triassic or Jurassic in age, as well as ophiolitic rocks that have been dated as Jurassic. In this study, we use the latter unit subdivision with a basal Veleta unit and overlying Mulhacén succession since it better reflects the initial protolith ages ( Figure 1). We utilize the previous chronostratigraphic framework established by Poulaki & Stockli (2022) and integrate it with microstructural observations and zircon and apatite U-Pb data to investigate the style and timing of deformation during subduction and underplating of the NFC units.
Previous work has already indicated a potential for tectonic relationships within the NFC. The first studies in the region described the contacts between these units as thrust faults (Egeler & Simon, 1969; García-Dueñas  Kampschuur & Rondeel, 1975), but other studies suggest the contacts could be mylonitic shear zones developed during nappe tectonics with thrust contacts in between the units that later cut through these contacts as low angle normal fault shear zones (e.g., Booth-Rea et al., 2005;Comas et al., 1999;Crespo-Blanc, 1995;Martínez-Martínez & Azañón, 1997;Martínez-Martínez et al., 2002, 2010. More recently, Sanz de Galdeano and Santamaría-López (2019) proposed that the Veleta-Mulhacén contact is stratigraphic and transitional in nature and argued for the lack of structural discrepancies between the units. The revaluation of the MDA calculations from Poulaki & Stockli (2022) indicated at least eight structural repetitions, and here we compared these ages with petrological data, temperature conditions, textural fabrics, and field observations across the NFC to reevaluate the potential for significant internal structural features related to subduction deformation and underplating processes ( Figure S5 in Supporting Information S1).

Cenozoic Tectono-Metamorphic Evolution of the Betic Cordillera
The tectonic processes that led to the closure of the Neo-Tethys ocean in the western Mediterranean and the Betic-Rif orogen are currently debated, with scenarios including (a) a northwest-dipping subduction zone (e.g., Bezada et al., 2013;Booth-Rea et al., 2007;Brun & Faccenna, 2008;Carminati et al., 2012;Chertova et al., 2014;Faccenna et al., 2004;Rosenbaum et al., 2002), (b) An east-dipping subduction zone followed by reversal and northwest-dipping subduction (Frizon de Lamotte et al., 2000;Rehault et al., 1984), (c) a south-east dipping subduction zone (Behr & Platt, 2012) or (d) contemporaneous southward and northward two-sided subduction Vergés & Fernàndez, 2012). Regardless of the initial direction of subduction, many studies agree on the subsequent rotation of the system with east-dipping subduction and westward trench retreat. Recent geophysical studies have imaged the Iberian Moho at depths of 65 km, which implies substantial crustal thickening (de Lis Mancilla et al., 2015). Seismicity data show a steeply dipping slab with earthquakes continuing to 120-150 km depth (e.g., Civiero et al., 2020;Heit et al., 2017). Geochemical and geophysical evidence of a Miocene volcanic arc related to subduction are evidenced in southeastern Spain (e.g., Booth-Rea et al., 2018;Casalini et al., 2022;Gómez de la Peña et al., 2020). Additionally, rocks with HP metamorphic signatures are exposed in the NFC, and the subcontinental Ronda peridotite is preserved in southern Iberia emplaced between crustal units of the Alpujarride nappe stack. Collectively, these lines of evidence indisputably confirm the involvement of a subduction zone in late Cenozoic convergence between Africa and Iberia; although, the timing of subduction initiation and evolution of the system remains obscured. During the Middle Miocene, westward trench retreat of the Betic subduction zone involved slab rollback and lateral migration to its current configuration (e.g., Moragues et al., 2021).

Methods and Results
In this study, we focus on the Cenozoic metamorphic and structural evolution of the NFC. We use a densely sampled data set across all structural levels of the NFC throughout the Betic Cordillera. Detrital zircon analyses of these samples were previously presented in Poulaki and Stockli (2022), which used provenance source signatures to determine the Paleozoic to Mesozoic depositional history and presubduction stratigraphy of the NFC. Here, we utilize this new stratigraphic framework and the spatial age constraints coupled with new apatite and zircon geo/thermochronology and microstructural analyses to study when and how units of the NFC were modified during subduction and exhumation.
Apatite petrochronology is a great tool to set relative temperature constraints on rocks that have experienced subduction-related metamorphism since its closure temperature of Pb has been determined by laboratory experiments to ∼360-550°C (Cherniak et al., 1991;Chew & Spikings, 2021;Smye et al., 2018;Watson et al., 1985). Additionally, apatite records deformation and fluid interactions and can be partially recrystallized and record multiple deformation and alteration events (e.g., Odlum & Stockli, 2020;Odlum et al., 2022). Even though few studies have investigated apatite U-Pb systematics in subduction rocks, Henrichs et al. (2018) showed that apatite is fully reset within its partial retention zone. Recent advances in zircon U-Pb geochronology have shown that zircon can decipher metamorphic events during the different stages of subduction metamorphism from prograde to peak to retrograde metamorphism by investigating metamorphic overgrowths (e.g., Kohn & Kelly, 2018;Poulaki et al., 2021;Rubatto, 2002). Lastly, quartz and feldspar microstructures with EBSD analyses complement these chronometers and provide relative temperature constraints.

Metamorphic Zircon
Of the 71 samples previously analyzed from the NFC in Poulaki & Stockli (2022), 33 of them yielded Cenozoic zircon overgrowths and are presented and discussed in this study ( Figure 1). These zircon grains were analyzed by laser-ablation inductively coupled plasma mass-spectrometry depth-profiling following the procedures of Marsh & Stockli (2015), which resulted in a continuous radiometric sequence from thin rims to zircon cores. We used GJ1 as the primary zircon standard (601.7 ± 1.3 Ma; Jackson et al., 2004) and Plešovice (337.1 ± 0.4 Ma; Sláma et al., 2008) as the secondary standard. Data reduction was carried out with IgorPro-based Iolite 3.4 software (Paton et al., 2010). By using depth-profiling techniques, we are able to manually differentiate and interpret age plateaus in the U-Pb concentrations and distinguish zones between cores, mixing ages, and rims of the zircon grains. After depth profiling, several grains were polished, and we collected Cathodoluminescence images (CL). Three grains were chosen for the collection of 2D zircon elemental mapping techniques as further described in Poulaki et al. (2021). We discarded zircon rim ages with more than 30% discordance between 206 Pb/ 238 U and 207 Pb/ 235 U ages and reported our data with 2σ propagated errors.

Zircon U-Pb
Calculated U-Pb ages of zircon metamorphic overgrowths range from 100 Ma to 12 Ma with the population majority spanning between 60 Ma and 17 Ma (Figures 2a-2c). From Sierra Nevada, 11 samples from the upper structural . Kernel Density estimate plots use the corrected ages by using the Stacey and Kramers (1975) common lead correction. "N" represents the sample number and "n" is the number of grains. Panels (a-c) have been grouped based on the location where the samples were collected. levels exhibited zircon metamorphic overgrowths. These samples are within the Tahal (Permian-Triassic protolith), Bédar-Macael (Carboniferous-Jurassic protolith), and uppermost Veleta (Devonian-Carboniferous protolith). Two main age peaks are found in the zircon age distributions clustered around Paleocene/Eocene ( Figure 2ai) and Miocene ( Figure 2aii). One sample, a quartz mica phyllite (18SSN13) from the uppermost Veleta, shows a single rim age of the Eocene. The absence of zircon rims is prominent in the central parts of the orogen and mostly within the Veleta unit as well as from four samples from the Mulhacén succession ( Figure 1).
From Sierra de los Filabres, 14 samples from the NFC exhibit a wider range of zircon rim ages with a minor age mode at ∼70 Ma. Clustering the Paleocene/Eocene ages together yields an average age of ∼55 Ma (Figure 2biii). The Miocene peak for these zircon metamorphic overgrowths is well defined with an average concordant age of ∼18 Ma ( Figure 2biv). The majority of the samples with zircon rims derive from the eastern Sierra de los Filabres exposes the uppermost units. Out of 14 samples that preserve Cenozoic zircon overgrowths, three of them have only Eocene rims (19SSF08,19SSF06,and 18SSF10) and two have only Miocene rims (19SSF12,19SSF07). The lack of Cenozoic zircon rims is observed in the largest part of the Veleta unit as well as in five samples from the Mulhacén succession. Similar to the samples collected from Sierra Nevada and Sierra de los Filabres, eight samples with Cenozoic zircon overgrowths from Sierra Alhamilla show two distinct peaks in the Eocene and Miocene. Most samples from the upper structural sections within the Bédar-Macael and uppermost Veleta (19SSA08) exhibit both generations of rims. Importantly, samples from the majority of the Veleta unit do not preserve any Cenozoic metamorphic zircon rims.
The spread of zircon metamorphic rim ages does not necessarily reflect the extent of the metamorphic event since they could also be influenced by discordance, inheritance, and common lead because of the thin nature of the overgrowths. Metamorphic rims older than 80 Ma are not further discussed since these data are sparse and do not represent a statistically robust portion of the collective zircon rim population. Furthermore, convergence between Africa and Iberia begins around ∼60-100 Ma; hence, it is unlikely that they are related to Cenozoic subduction metamorphism, which is the focus of our study.

Zircon CL Images and 2D U-Pb Maps
To better constrain the spatial architecture of zircon overgrowths, we collected CL data to image zonation patterns and 2D LA-ICP-MS elemental maps to constrain their radiometric ages. Several representative grains from three orthogneisses (19SSA10,19SSN12A,18SSF08) with Eocene and Miocene rims revealed by depth-profiling data were selected for CL analyses. These grains were hand-picked, rotated, and mounted to double-sided sticky tape, covered with epoxy, and polished for analysis. CL images reveal that Eocene rims are generally thin with homogeneous textures (Figure 3a and Figure S1 in Supporting Information S1). These rims are darker than their corresponding magmatic cores and follow the oscillatory zoning geometry of the grain cores but preserve different textural morphology and are overall euhedral. In contrast, Miocene rims appear thicker and more heterogeneous with complex porous spongy textures ( Figure 3b and Figure S1 in Supporting Information S1). These rims are asymmetric with variable thickness around the grain edges and are brighter than their corresponding cores.
Three representative grains from sample 19SSA10 were selected for 2D elemental mapping. These analyses provide a 2D map of age constraints throughout the zircon grain and thus reveal their isotopic morphology. In Figure 3a (19SSA10_16) the Eocene rim is only ∼10 μm thick and has low (<0.05) Th/U with high U concentrations. The core of this grain is Permian and clear zoning is observed in the corresponding CL image with low U (∼530 ppm) concentration and high Th/U. Zircon grains 19SSA10_52 and 19SSA10_60 (Figures 3b and 3c) have Miocene rims enclosing a Permian core. These rims are as thick as 50 μm in places with lower Th/U (∼0.01) relative to their cores (>0.06). CL images show that these rims appear to have variable textures, with some regions growing along the oscillatory zoning, but others seemingly intruding and overprinting core zonation with porous structures.

Apatite U-Pb Geochronology and CL Imaging
Apatite U-Pb data were collected from samples from the lower-and upper-unit successions and analyzed both in situ on thin sections as well as on grain separates. Prior to in situ analyses, we collected element dispersive Elemental Dispersive Spectroscopy (EDS) maps to measure Ca, P, Si, and Ti and identify apatite within our thin sections. Afterward, we aligned the EDS map to the laser stage to target the mapped apatite grains. A 30 μm spot with a total of 200 ablation shots was used for in situ thin section analyses, while a 40 μm spot with 300 ablation shots was used for apatite grain separates. For selected grains, we performed 2D elemental maps using 15 μm and following the same methodology as for the zircon maps. For all three methods, MAD apatite (Thomson  Tera and Wasserburg, 1972;Vermeesch, 2018).
In our samples from the NFC, we distinguish four distinct groups of apatite grains based on their ages and textures. The first group is defined by late Paleozoic and early Mesozoic U-Pb ages, which are older than subduction initiation and convergence between Africa and Iberia. Samples in this group are solely from the Veleta unit, with two in situ analyses (19SSN01,19SSN02) and two from grain separates (18SSN09, 18SSA02). Graphitic mica-schist samples (18SSN09, n = 13; 18SSA02, n = 44) with Carboniferous MDAs yielded apatite U-Pb ages of 379.1 ± 11.4 Ma and 307.6 ± 2.6 Ma, respectively, and a common Pb composition of ∼0.84 (Figures 4a and 4c). A Carboniferous graphitic quartz mica schist (19SSN01, n = 44) yielded in situ apatite ages of 152.3 ± 7.0 Ma. Similarly, in situ apatite from a metabasite within the Veleta unit yielded a Jurassic age of 197.0 ± 6.0 Ma (Figure 4b;   Poulaki & Stockli, 2022). CL images from these apatite grains show an overall homogeneous grain structure with minor patches of bright spots on grain surfaces ( Figure 4 and Figure S2 in Supporting Information S1).
The second group is characterized by Eocene apatite cores and rims. This group contains two Permian tourmaline-bearing orthogneisses from the Mulhacén succession in the Sierra Nevada (19SSN12A, n = 484; 19SSN20, n = 88) and were previously found to have zircon crystallization ages of ∼275 Ma (Poulaki & Stockli, 2022). Apatite U-Pb ages from grain separate from these samples are Eocene with ages at 45.0 ± 2.2 and 42.5 ± 5.8 Ma, respectively, and a common Pb composition of ∼0.7 (Figures 5a and 5b). Due to the complicated structure revealed by the CL images and the high common Pb composition in sample 19SSN12A, we exported and plotted the data in 3-s increments as described in Odlum & Stockli (2020). CL images from these grains exhibit bright irregular cores with darker rims. Additionally, a network of fractures crosscutting in irregular orientations is observed on many of these grains (Figure 3d). Hence, it appears that these bright fractures postdate apatite rim formation. With our depth-profiling method for these grains, we can subdivide the raw data acquisition window into distinct plateaus, which yields an Eocene age for the apatite core. However, we were not able to obtain ages for the rims due to low U and high concentration of common Pb. In our 2D U-Pb maps, the apatite core and rims can be easily differentiated by differences in the U and Th concentrations (Figures 3d and 3e). In both grains from samples 19SSN12A and 18SSN20, the cores are highly enriched in U (∼60-100 ppm) in comparison to the rims (∼10-30 ppm). In contrast, Th concentration in grains from sample 19SSN12A is the same in rims and cores (∼20 ppm) but higher in cores (∼20 ppm) than rims (∼5 ppm) in grains from sample 19SSN20 (Figures 2d and 2e).
The third group of apatite grains were collected from metasedimentary rocks from the Mulhacén succession at Sierra Nevada and are distinguished by Miocene apatite U-Pb ages. A Triassic protolith garnet quartz mica schist (18SSN14, n = 123) was analyzed in situ and yielded an age of 9.4 ± 4.9 Ma and a common Pb composition of 0.8. A Devonian quartz rich schist (19SSN11, n = 341) yielded a U-Pb apatite age of 16.2 ± 3.7 Ma and a common Pb composition of 0.83. Similar to sample 19SSN12A, we exported these data in 3-s increments. CL images from these grains show an amorphous core with many dark and light patches and bright rims around the cores (Figures 5c and 5d).
The fourth and final apatite group is distinguished by partially reset apatite. Grains for this group are from samples collected from Sierra Alhamilla and Sierra de los Filabres. Our data show that these grains are partially recrystallized, evidenced by large dispersion on individual analyses within the Tera-Wasserburg plots (Figures 5e and 5f). By isolating the different ellipse clusters within the age spectra, we can reconstruct the oldest and youngest events the apatite has recorded from the grouped lower intercepts. A calcite mica quartz schist (19SSF09) from Sierra de los Filabres displays a lower youngest intercept of 47.1 ± 16.7 Ma and a lower oldest intercept of 270.1 ± 11.1 Ma (Figure 5e). Orthogneiss samples from Sierra Alhamilla (18SSA10) yielded a younger lower intercept of 48.6 ± 2.2 Ma and an older lower intercept of 293.2 ± 12.6 Ma ( Figure 5f). CL images from these grains show a darker core with lighter patches and a very bright rim with a porous, spongy texture ( Figure 5f). The thickness of rims varies from 1 to 2 μm up to 20 μm. Last, a Devonian quartzite (19SSN07) in the southern Sierra Nevada exhibited an apatite age of 40.9 ± 17.1 Ma ( Figure S3 in Supporting Information S1). Due to the large error attributed to low U concentration and high amounts of common Pb, we do not lump this sample into the four apatite groups.

Microstructural Analyses
Since the relationship between temperature and deformational patterns expressed in quartz and feldspar microstructures is well established, we conducted detailed observations on samples from Sierra Nevada to complement apatite U-Pb thermo-chronometer data. Collectively, these tools provide estimates of the relative temperatures (∼400-550°C) that these rocks experienced during their polymetamorphic evolution and place timing constraints on strain accommodation during these events. Microstructures were qualitatively assessed with petrographic microscopes and quantified using EBSD techniques. All samples for these analyses were cut perpendicular to the foliation and parallel to the lineation. EBSD mapping was performed at the Geomaterials Characterization and Imaging Facility at the University of Texas at Austin. We used AzTec software to clean the raw quartz and feldspar EBSD data by removing wild spikes and plotted the data using the MTEX software package (Bachmann et al., 2010). In this section, we describe the textural observations in structural order from the lower to upper units of the NFC.
Two graphitic quartz mica schists from the Veleta unit (19SSN01,19SSN09) show similar microstructural characteristics. Both exhibit ∼1 mm euhedral quartz crystals within microlithons (Figures 6a, 6b, 6e, and 6f). Misorientation maps from these samples record limited internal strain, with a few quartz grains preserving some subgrains and minor undulatory extinction and minor evidence of grain boundary migration. However, the majority of the quartz crystals in these samples show no evidence of intracrystalline plasticity. We also found evidence for normal grain growth with 120° angles in quartz grains. Mica in these samples are aligned along the quartz ribbons. EBSD data show that these samples have absent or weak Lattice Preferred Orientation (LPO) as defined by their Misorientation index (M) (M = 0.0038 and M = 0.0068; Figures 6c, 6d, 6g, and 6h), but there is a Shape Preferred Orientation (SPO) parallel to the foliation. In sample 19SSN01, albite crystals contain graphitic inclusions that define microfolds ( Figure S4 in Supporting Information S1) and show some subgrains and minor undulatory extinction but no evidence of major plastic deformation within the grains. Albite grains contain inclusion graphite trails and exhibit evidence of rotated textures relative to the main foliation orientation.
Samples from the higher parts of the Veleta unit contain garnet porphyroblasts, in contrast to the lower sections that do not contain garnets, as documented by this work and previous studies (e.g., Behr & Platt, 2012;. A quartz mica schist with feldspar and minor garnet and a quartz-rich schist (19SSN03) collected from within the uppermost Veleta close to the Veleta-Tahal contact have quartz crystals with distinct layers of elongated ribbons with grain boundary migration and bulging (Figures 7a and 7b). Additionally, quartz forms large ameboid grains and has a strong LPO (M = 0.1036), as shown by the mis2mean maps (Figures 7c and 7d), which are indicative of moderate temperatures (∼500°C; Passhier & Trouw, 2005). Similar structures have been observed by Behr and Platt (2013) in Sierra Alhamilla. Albite porphyroblasts in this sample are euhedral with minor subgrains but overall do not show evidence of recrystallization (Figure 7b). A Devonian quartzite from within the Bédar-Macael unit (19SSN11) also exhibits a moderate LPO (M = 0.0833) with bulging in quartz and grain boundary rotation (Figures 7e, 7f, and 7g). However, the presence of feldspar grains in between quartz grains in this sample could potentially hinder the development of quartz LPO fabrics.
Within the Tahal schist, a garnet mica schist (18SSN14) preserves significantly different microstructures than previous samples from lower stratigraphic domains. Garnet porphyroblast growth predates the current foliation, and they have well-defined inclusion trails with quartz, mica, and minor rutile. Quartz within the main foliation shows evidence of grain boundary migration and recrystallization. Albite crystals show significant dynamic recrystallization and breakdown due to mica growth and subgrain rotation (Figures 7h and 7i), with evidence of subgrains indicating temperatures >550°C. Within the uppermost Bédar-Macael unit, two tourmaline-bearing orthogneisses were further investigated (19SSN12A, 18SSN20, 18SSF08). Feldspar grains have discontinuous undulatory extinction, subgrains, and bulging, but no dynamic recrystallization is observed (Figure 6j). Quartz within these samples shows grain boundary migration with ameboid structures and many subgrains (Figure 6j).

Stratigraphic Rearrangement: Current Structural Position of NFC Units
Subduction complex rocks are highly deformed and often have experienced multiple metamorphic episodes, causing much of their initial crustal anatomy and stratigraphic architecture to be overprinted and obscured. Some rocks from subduction complexes appear to have been deformed and/or underplated in a mélange style with significant levels of lateral flow and mixing, resulting in a complex loss of presubduction context (e.g., Angiboust et al., 2013;Cloos & Shreve, 1988). These HP rocks are often accreted to the overriding plate in a complex fashion and subsequently exhumed, involving orogenesis, uplift, erosion, and syn-and postorogenic extension. Hence, constraining the temporal and spatial scales of progressive deformation and underplating as well as metamorphism along the deep subduction interface is challenging solely based on large-scale field observations. Recent breakthroughs combining geochronology with detailed field observations in the Cyclades of the Hellenic subduction zone revealed large-scale stratigraphic coherency of these deformed rocks with structural imbrication and stacking of tectonic slivers (Kotowski et al., 2022;Poulaki et al., 2019;Seman, 2016). In this study, we leverage this approach by combining detailed geochronologic analyses and utilizing the presubduction stratigraphic reconstructions of the NFC based on MDAs established by Poulaki & Stockli (2022). Their work showed that the provenance signature of the NFC is remarkably stable throughout its depositional history and largely records unroofing of late Paleozoic Variscan basement and early Mesozoic Alpine Tethyan rifting. The consistent younging of MDAs within tectonic slivers, regional lithological correlations, and dating of igneous orthogneisses provide strong evidence that these zircon ages provide a robust chronostratigraphic framework for the NFC. For this study, we interpret MDA patterns from 12 regional transects across the NFC to assess how these rocks have been structurally rearranged and stacked during Cenozoic subduction and exhumation. Throughout the structural column of the NFC, we identified at least three MDA reversals indicative of older on top of younger rocks, and eight total reversals within our transects (Figure 8). These reversals show that the original stratigraphic architecture of the units has been modified, requiring a structural explanation, such as the presence of either recumbent folds or thrust fault imbrication. While small-scale folds are present, large-scale recumbent folds (Arend Zevenhuizen, 1989;García-Dueñas et al., 1988;Martínez-Martínez et al., 2002, 2010 should result in age reversals across the fold axis and symmetry of lithological and age progressions from old to young and to old again. Instead, we consistently observe younging MDA estimates within each separated sliver and abrupt changes back to old ages at lithological boundaries. These relationships are best explained by the presence of discrete thrust faults that imbricate coherent and upright tectonic slivers, forming a nappe stack and preserving original stratigraphic packages that are young within each individual fault-bounded sliver. Additionally, we find consistent foliation patterns across the MDA-defined sliver boundaries (Figure 8; Ruiz-Fuentes & Aerden, 2018) but differences in metamorphic pressure/temperature conditions and metamorphic and maximum depositional ages (this study; e.g., Santamaría-López et al., 2019). Collectively, these lines of evidence argue for the presence of multiple thrust repetitions and stacking within the NFC. Below, we integrate the MDAs, microstructures, field observations, and metamorphic ages to discuss where, when, and how these large-scale (∼100s m) deformation features formed and the implications for subduction underplating processes.
The majority of the Veleta unit exhibits monotonous MDAs ranging from Early to Late Carboniferous and we interpret this unit as a thick coherent Variscan basement block that did not experience any significant post-Variscan internal imbrication and hence largely preserves its presubduction architecture. However, along the uppermost parts of Veleta unit in the Sierra Nevada and Sierra de Los Filabres, we find Cambrian/Devonian rocks intercalated with Carboniferous strata, and in some cases, sandwiched between Permian strata (Figures 8a, 8b, and 8e). This age reversal corresponds to the base of the Calar-Alto unit mapped by Martínez-Martínez et al. (2010). The observed Devonian on top of Carboniferous relationships could be explained as thrusts or large folds formed during the Variscan orogeny. However, we observe similar relationships with Devonian atop Permian rocks in the western Sierra Nevada and southern Sierra de los Filabres. Additionally, the first age reversal above the Veleta coincides with the position of the Dos Picos shear zone, where hotter rocks where emplaced over lower-temperature ones at a late stage in the metamorphic evolution of the NFC complex (Augier, Agard, et al., 2005;B. Li & Massonne, 2018). Together, these evidence show that these tectonic contacts postdate the Variscan deformation in the Carboniferous and imply thrust imbrication during Cenozoic convergence and subduction (Figures 1 and 8).
Moving up section, the relationship between the overlying Mulhacén succession and the Veleta unit is spatially variable. While in most locations, we observe a continuous depositional relationship from Carboniferous to Permian rocks (Figures 1 and 8f), supporting similar findings by Sanz de Galdeano and Santamaría López (2019), along the western domain of Sierra Nevada and Sierra Alhamilla, Carboniferous rocks of Veleta are in direct contact with Triassic strata of the Tahal unit of the Mulhacén Succession (Figures 1 and 8a). This hiatus could be interpreted as a synrift unconformity in the depositional record, although a more complicated structural relationship between the Veleta and Mulhacén units has been proposed (Ruiz-Fuentes & Aerden, 2018). Alternatively, this age relationship could also likely be related to extensional shearing during late-stage exhumation .
The Mulhacén succession exhibits the most variability in lithologies and MDAs, with mixed metabasites, marbles, schists, and orthogneisses, and ages ranging from Early Permian to Early Jurassic. There is ample evidence for older on younger MDA relationships within the Mulhacén succession, indicating thrust relationships and pervasive structural imbrication. Several locations show slivers of Permian rocks wedged above Triassic but below Devonian strata (Figures 8a and 8b). Similar tectonic slivers have been previously identified in the eastern Sierra de los Filabres and were interpreted as three distinct thrust sheets with Paleozoic on top of Mesozoic cover (de Jong, 1993a(de Jong, , 1993bPorkoláb et al., 2022). Our new data confirm this interpretation and show Permian rocks of the lower unit thrust above Triassic strata (Tahal formation; Figure 8d). We also observed similar imbrication patterns in the upper structural portions of the Sierra Nevada. Comparable age reversal patterns are present in the western and southwestern Sierra Nevada, where Permian rocks are structurally juxtaposed against Triassic and Devonian slivers. The dominant S2 foliation in all slivers consistently dips direction toward N/NW in West Sierra Nevada and NNE in Sierra de los Filabres (Figures 8d and 8e). Ruiz-Fuentes and Aerden (2018) showed that S1 and S2 foliations associated with subduction fabrics are slightly more variable within the Mulhacén and along the Mulhacén/Veleta contact, whereas they are more consistent throughout the Veleta unit, supporting our observations that Mulhacén has more internal structural complexity. Thrust relationships, similar to those in the NFC, have also been identified in the adjacent Alpujárride Complex by field observations (Balanyá et al., 1997;Booth-Rea et al., 2005;Rossetti et al., 2005) and recent geochronologic data (Poulaki & Stockli, 2022), showing Carboniferous graphitic schist structurally overlaying Permian phyllites and schists. Hence, both the NFC and Alpujárride units have experienced subduction imbrication and subsequent underplating and exhumation.

Pre-Cenozoic Metamorphism of the Veleta Basement
The integration of zircon and apatite U-Pb and microstructural analyses strongly suggests that the NFC Veleta unit primarily records presubduction deformation and metamorphism during the late Paleozoic Variscan orogeny, despite being subducted and exhumed during the Cenozoic. In general, apatite grains are rare in the lithologies of the Veleta unit, and only four samples yielded apatite grains. Two samples from Veleta in Sierra Nevada and Sierra Alhamilla gave Carboniferous apatite U-Pb ages (Figure 4). While these metasedimentary protoliths are characterized by Carboniferous MDAs and cosmopolitan detrital zircon spectrum spanning 3 Ga to 300 Ma (Poulaki & Stockli, 2022), the Carboniferous/Devonian apatite U-Pb ages appear to record Variscan tectonism and metamorphism.
A metabasite and graphitic mica schist from the Veleta unit in the southern Sierra Nevada both yielded the earliest Jurassic apatite U-Pb ages (Figures 4a and 4b). The metabasite age likely represents the timing of magmatic crystallization and has been interpreted as recording dike emplacement as part of the Central Atlantic Magmatic Province and coeval with early Jurassic rifting (Poulaki & Stockli, 2022). The Late Jurassic apatite U-Pb age from the Carboniferous country rock records subsequent rifting and break-up during the formation of the Alpine Tethys. This age is consistent with previously reported zircon ages from metabasites within the NFC (Puga et al., 2002) and has been attributed to Jurassic rifting in the Western Mediterranean following the Variscan Orogeny (e.g., Saspiturry et al., 2019). Interestingly, Cenozoic metamorphic zircon rims in the Veleta unit are sparse or entirely absent (Figures 1, 8c, and 8f).
In terms of metamorphic paragenesis, the lower part of the Veleta unit is dominantly composed of quartz, feldspar, muscovite, biotite, chlorite, and rutile, consistent with these rocks not having exceeded midgreenschist facies metamorphic conditions. Despite evidence of later static recrystallization with some polygonal quartz grains, EBSD data show a very weak or absent LPO fabric in quartz grains. While experimental work has shown that quartz LPO is commonly preserved even when rocks undergo static recrystallization (Heilbronner & Tullis, 2002), the evidence for some static recrystallization of quartz does not support the annealing or complete removal of a preexisting fabric. Additionally, Veleta samples are characterized by the lack of evidence for intracrystalline plasticity in quartz and recrystallization of feldspar (Figures 6a-6h). Similar microstructures have been observed in the Veleta unit in the Sierra de los Filabres (González-Casado et al., 1995) and Sierra Alhamilla . The latter study attributed the lack of LPO and presence of SPO to pressure-solution creep during subduction and exhumation with the presence of free water at estimated temperatures of 490-530°C (Behr & Platt, 2012).
Our new apatite U-Pb data from the majority of Veleta yield robust preCenozoic ages, suggesting that this unit did not experience temperatures >450°C during subduction-related metamorphism and escaped major fluid alteration necessary for apatite recrystallization. We propose that the Veleta unit experienced high-temperature metamorphism (>∼500°C) during the late Paleozoic Variscan tectonic events but remained largely unaltered during Cenozoic subduction. Presubduction metamorphism of the Veleta unit may have affected the rheology and deformation behavior of this section of the NFC during Cenozoic subduction compared to the overlying post-Variscan NFC strata, which show evidence for younger apatite ages and hence higher temperature subduction metamorphism.

Cenozoic Metamorphism During Subduction and Underplating
Systematic LA-ICP-MS depth profiling analyses revealed two primary metamorphic events in the Eocene and Miocene recorded by metamorphic zircon rim ages and metamorphic apatite grains (Figures 2 and 5). A Miocene age for metamorphism of the NFC has been widely reported in numerous studies; however, Eocene metamorphism remains controversial. This marks the first study to propose a polyphase metamorphic evolution for the NFC from multiple geochronometers and consistently resolve the multistage Cenozoic evolution of the region.

Eocene Metamorphism
An Eocene HP/Low-Temperature (HP/LT) metamorphic phase in the NFC has been suggested based on 40 Ar/ 39 Ar analyses (Augier, Agard, et al., 2005;Monié et al., 1991;Porkoláb et al., 2022), while other studies have disregarded these data and attributed them to excess Ar because they do not fit evidence for younger Miocene metamorphic ages interpreted as a HP/LT event (Behr & Platt, 2012, de Jong, 2003De Jong et al., 2001;Kirchner et al., 2016;Platt et al., 2006). However, B. Li and Massonne (2018) identified Eocene monazite ages in the Mulhacén succession of the Sierra Nevada that support the Eocene 40 Ar/ 39 Ar data and suggest that the subduction of the NFC must have occurred prior to the Miocene metamorphic event. This was corroborated by Aerden et al. (2022), who reported a late Eocene garnet Sm-Nd age (∼35 Ma).
This study presents robust apatite and zircon overgrowth U-Pb ages that document a well-defined Eocene/Paleocene metamorphic event between ∼60 and 35 Ma (Figures 2, 5a, 5b, 5e, and 5f). On the basis of these data, we suggest that this event represents a metamorphic event with temperatures that exceeded ∼450°C in the Betic subduction zone. The Eocene metamorphic event is also recorded by apatite groups 2 and 4 ( Figure 2). Zircon rims are observed at the upper parts of the structural column within the Mulhacén succession and the uppermost Veleta unit. In most cases, zircon rims are euhedral while still following the oscillatory zoning of their cores ( Figure 3a) and have Th/U ratios consistently less than 0.1, which argues for their metamorphic nature ( Figure  S6 in Supporting Information S1; e.g., Rubatto, 2002;Williams, 2001). With our depth-profiling techniques, we were able to analyze the very thin rims (∼5 μm) and recover a metamorphic population that previous studies overlooked. Eocene rims are not observed throughout most of the Veleta unit and, in some cases, on samples from the hanging wall of the thrust faults (Figure 1), which confirms their subduction-related genesis and supports differential metamorphic conditions across imbrication thrusts and argues for a post-Eocene age of fault activity.
Eocene apatite grains from Bédar-Macael orthogneisses (Group 2) exhibit bright cores and dark rims (Figures 5a  and 5b). Our apatite U-Pb depth-profiling data and the 2D apatite U-Pb mapping indicate that these bright core zones correspond to high U concentrations (Figures 3d and 3e). Previous studies have shown that recrystallized apatite rims have significantly lower U concentrations than their cores (Henrichs et al., 2018;O'Sullivan et al., 2020), which is in agreement with our 2D maps. Given the high U cores, we are confident that the Eocene ages correspond to the core of the grains, and hence the original magmatic apatite cores must have been completely reset. Due to the low U concentration and high common Pb composition of the rims, we are unable to obtain a robust U-Pb age, but their formation must postdate the Eocene metamorphism. These samples also contain quartz ameboid structures with grain boundary migration and feldspar crystals with evidence of bulging. Together, the microstructural observations and apatite ages indicate that these rocks exceeded ∼450°C during subduction metamorphism, which agrees with previously published P/T estimates (Augier, Agard, et al., 2005;Booth-Rea et al., 2015;Ruiz-Cruz et al., 2015).
Apatite grains in Group 4 from Sierra de los Filabres and Sierra Alhamilla show a large dispersion of ages in the Tera-Wasserburg space with lower-intercept ages ranging from Carboniferous to Eocene (Figures 5e and 5f). We interpret this age dispersion as resulting from partially recrystallized apatite, where some zones of the crystal lattice preserve the Eocene metamorphic event. Similar patterns of partially recrystallized apatite have been observed in mylonitic zones in the Pyrenees (Odlum & Stockli, 2020;Odlum et al., 2022), and the youngest ages have been thought to record the most recent deformation event. Thin sections from our samples show evidence of fluid-rock interactions with feldspar replacement by sericite (Figure 6i). Because apatite grains from Sierra Nevada are fully reset, the partially recrystallized apatite from Sierra de los Filabres and Sierra Alhamilla might indicate that these rocks experienced lower peak temperatures. Similar temperature variations between the Betic orogens have been previously observed in both the NFC and Alpujárride with higher temperatures in the west and lower in the eastern domains (Bessière et al., 2022;Platt et al., 2006). Additionally, different P/T conditions have also been observed among units in the same localities, which results in a discrepancy among P/T work done in the area as compiled by Santamaría-López et al. (2019).
Our new data unequivocally show strong evidence for an Eocene metamorphic event and are coeval with Eocene metamorphic ages found by previous studies. Collectively, the textures we observe in apatite and zircon rims, the closure temperature of apatite at ∼450°C, and the general absence of zircon rims in the structurally deepest Veleta unit and across the thrust contacts confirm that metamorphism is subduction-related and not formed prior to underthrusting.  performed a detailed P-T study from syn-kinematic whitemica ages that conclusively showed that HP conditions existed throughout the Eocene. Therefore, it is most likely that the growth of our zircon rims and resetting of apatite grains can be attributed to HP metamorphic conditions during the Eocene, although additional geochemical and in situ data in HP assemblages are needed to holistically evaluate the nature of Eocene metamorphism.

Miocene Metamorphism
The Miocene metamorphic event has been well documented in the NFC and has been generally interpreted as an HP/LT metamorphic event with age constraints spanning from ∼20 to 13 Ma (Gómez-Pugnaire et al., 2012;Lopez Sanchez-Viscaino et al., 2001;Platt et al., 2006). We find a distinct Miocene zircon rim population spanning ∼30-10 Ma, in good agreement with ages from previous studies. The Miocene zircon metamorphic overgrowths show a different morphology than the Eocene ones, with porous structures and thick, uneven rims around the cores of the zircon that suggest a metasomatic nature for these rims (Figures 3b and 3c) (e.g., Poulaki et al., 2021). The majority of samples with Miocene zircon rims also have Eocene zircon rims and are found only in the upper structural levels of the NFC. The two samples that yielded Eocene core apatite ages (Group 2) also yielded Miocene zircon rims (Figures 5a and 5b). Hence, it is possible that the Group 2 apatite rims from these grains are Miocene in age and coeval with zircon rim formation, but we were not able to obtain an age due to the high common lead composition and low U concentrations (Figures 3d and 3e).
Additionally, we found Miocene (∼15 Ma) apatite ages (Group 3) in the Mulhacén succession in the Sierra Nevada (Figures 5c and 5d). These grains have complicated internal textures with multiple mixed zones of bright rims and patchy cores (Figure 5d and Figure S2 in Supporting Information S1). Even though we depth-profiled these grains, the complicated internal structure does not allow us to clearly differentiate if the rim or core of the apatite corresponds to the Miocene age we obtain; nevertheless, the cores must be reset in the Miocene since we do not observe a dispersive age spectrum as seen in Group 4. Microstructures from this sample (19SSN11) show a strong maximum in the center, but a large opening angle, and feldspar does not show recrystallization, possibly indicating lower temperatures (Figures 7f and 7g). Another sample with Miocene apatite contains garnet, and the quartz grains are elongated and exhibit grain boundary migration and feldspar recrystallization (Figures 7h  and 7i). The recrystallization of feldspar with subgrains and subgrain rotation are indicators of temperatures at ∼550°C. However, the ages of apatite do not necessarily record the precise timing and magnitude of these temperatures since apatite could reset under lower temperatures (∼450°C) and be influenced by the presence of fluids that these rocks experienced during Early Miocene reheating and rapid exhumation. Hence, recrystallization of apatite might have occurred at lower temperatures and postdate the recrystallization of feldspar grains.
Metamorphism in the NFC during the Miocene is complex, and more recent studies have challenged the classic HP/LT interpretation of this event. Various scenarios include stages of late reheating (Bakker et al., 1989;Booth-Rea et al., 2015;B. Li & Massonne, 2018;Vissers et al., 1995), decompression (Augier, Agard, et al., 2005;Ruiz-Cruz et al., 2015), and/or fluid involvement that contributed to new mineral growth such as apatite, zircon, and garnet. This reheating stage might also be the reason for the confusion in the literature regarding the nature of Miocene metamorphism and whether it reflects an HP/LT or HT/LP event. Studies dating garnet argue for HP/LT metamorphism in the Miocene (e.g., Kirchner et al., 2016;Platt et al., 2006). If the Eocene metamorphic event also represents HP/LT metamorphism, that could possibly indicate diachronous metamorphism of the NFC units, with peak metamorphism experienced at different times. Our data show strong evidence for metasomatic apatite and zircon formation during the Late Miocene (∼15 Ma). These ages are slightly younger than the Early Miocene ages attributed to HP/LT metamorphism by previous studies (Kirchner et al., 2016;Platt et al., 2006). Therefore, we suggest that our Miocene chronometers may record both HP/LT peak metamorphic conditions in the Early Miocene followed by reheating and fluid-rich metamorphism coeval with exhumation in the Late Miocene. Previous work has identified at least five stages of fluid pulses starting in the Miocene, with the initial pulse corresponding to metamorphic fluids related to reheating and dehydration reactions (Dyja et al., 2016). The involvement of fluid pulses in the Late Miocene suggests that the NFC was structurally positioned above a dehydrating slab and actively exhumed in the overriding plate at this time. The spatial distribution of observed Miocene metamorphism only in the upper unit may indicate that fluid-rock interactions were preferentially localized along the upper units where the weak strata initially deformed by thrust faulting and later reactivated as extensional shear zones may have facilitated fluid flow (Martínez-Martínez et al., 2010;Porkoláb et al., 2022). The distribution of new metamorphic ages found in our study also aligned well with previous descriptions of the metamorphic history and structural mapping by Martínez-Martínez et al. (2010).

Reconciling Temperature Estimates and Chronological Constraints
The majority of the Veleta unit does not show any MDA duplications and has consistent foliation fabrics, arguing against widespread structural imbrication. Microstructural analyses show minor subgrains with no quartz LPO or feldspar deformation, indicating either static recrystallization or pressure solution mechanisms (Passchier & Trouw, 2005;Rutter, 1983;Stöckhert, 2002). Studies from the exhumed Arosa zone in the Alps have shown that pressure solution is dominant under temperatures less than ∼400°C (Condit et al., 2022). These microstructures, combined with pre-Cenozoic apatite U-Pb ages, suggest that the majority of Veleta experienced lower temperatures than the rest of the NFC, likely less than ∼450°C. In contrast, MDAs and petrological relationships indicate older-on-younger relationships at the upper structural parts of the NFC, including the Mulhacén succession and the uppermost Veleta unit. At these structurally higher levels, we observed Cenozoic metamorphic zircon rims and Cenozoic apatite grains (Figures 1, 8c, and 8f). In the upper section of the Veleta unit, a higher degree of metamorphism is evidenced by the appearance of garnets, intracrystalline plasticity in quartz and feldspar, and reset Cenozoic apatite ages. The Mulhacén succession likely recorded the highest overall temperatures among the NFC inferred by intracrystalline plasticity in both quartz and feldspar grains, with higher temperatures in Tahal exhibiting dynamic recrystallization in feldspar and lower temperatures in the Bédar-Macael with bulging in feldspar grains. In Sierra Nevada, these upper units exceeded ∼550°C with reset apatite, grain boundary migration in quartz, and recrystallization of feldspar grains.
Previous studies observed this apparent upward increase in the metamorphic record and attributed it to an inverted thermal gradient in the Betic subduction zone. Specifically, Behr and Platt (2013) argue that the upper parts of Veleta unit experienced deformation with pressure solution creep under relatively lower temperatures of 490-530°C (Behr & Platt, 2012), while the Mulhacén succession deformed at higher temperatures under dislocation creep locally up to 623 ± 68°C (Behr & Platt, 2012). Santamaría-López et al. (2019) demonstrated the similarity of P/T paths among the different NFC units with higher temperatures at the upper parts and also argued for an inverted thermal gradient rather than tectonic relationships between the units. Another scenario proposed that overheating occurred during this time due to contact from the overriding hot Alpujárride complex, causing heat and fluid transfer (Aerden et al., 2013;Bakker et al., 1989). We suggest that differences in the temperatures and deformation mechanisms are instead best explained by variations in the peak subduction depths throughout the NFC. This interpretation is supported by thermal grading from fully reset apatite in the West to partially reset apatite grains in the East within the same level of the Mulhacén Succession along the strike.
The Carboniferous and Jurassic metamorphic ages of the Veleta unit reveal that it has previously experienced metamorphism and hence dehydration prior to entering the subduction system. Additionally, the apparent unconformities with missing Permian strata and overall thicker stratigraphic sequence in the Veleta unit could indicate that it was deposited in a more proximal domain to the Variscan Orogeny and Pangea rifting events than the overlying NFC units (Poulaki & Stockli, 2022) despite subsequent minor thinning due to extension and detachment faulting during the latest Miocene. These conditions would result in a more buoyant succession that would provide resistance to deep subduction. Accordingly, previous work estimated pressures in the Veleta unit 2-8 kbar lower than the overlaying Mulhacén succession B. Li & Massonne, 2018). Our data argue that the lower Veleta basement could not have exceeded ∼450°C during subduction since apatite grains do not record any Cenozoic metamorphism. Collectively, we propose that the monotonous Veleta unit acted as a thick and homogeneous basement suite during the subduction of the NFC. The buoyant nature of Veleta may have provided a mechanism to trigger underplating of the NFC units and transfer to the overriding plate of the Iberian margin.
The similarity in provenance data between NFC and Alpujárride suggests that these units may have been laterally continuous throughout their depositional history but experienced different levels of Cenozoic metamorphism. These differences can be explained by the subduction of the NFC to greater depths beneath the Iberian margin than the Alpujárride Complex. In this study, the new age constraints from the NFC indicate that Alpujárride and the NFC were subducted at similar times. Since Alpujárride is currently on the hanging wall of low-angle normal faults, Alpujárride must have been subducted first, followed by the NFC in a lateral continuation of the same subduction system as recently described by several studies (Aerden et al., 2022;Porkoláb et al., 2022;Poulaki & Stockli, 2022). Previous models suggest that the subduction of the NFC occurs much later in the Miocene (e.g., Van Hinsbergen et al., 2020); however, our data, along with recently published data (Aerden et al., 2022;B. Li & Massonne, 2018;Porkoláb et al., 2022), suggests that the subduction of the NFC is likely to have commenced around the Early Eocene (Figure 9). New precise metamorphic ages in this study show that the Betic subduction zone initiated at least since the Eocene and progressively subducted and underplated units of the NFC. Our zircon rim age constraints on the timing of subduction zone metamorphism give an approximate timing of initial underthrusting around ∼60 Ma, which reconciles the controversy regarding subduction initiation and the large range of metamorphic ages spanning from Eocene to Miocene.

The Buoyancy of the Crystalline Basement (Veleta)
Geochronological data, combined with the MDAs along other transects, indicate that the Veleta unit did not form any large duplications or repetitions during subduction or exhumation and has also experienced metamorphism prior to subduction. We suggest that the Veleta unit behaves as a thick buoyant basement during the subduction of the NFC due to its large thickness and its prior dehydration from multiple previous metamorphic events during the Variscan orogeny and subsequent Jurassic rifting (Figure 10b). Prior metamorphism significantly contributes to the homogenization of the unit and may anneal some of the previous stratigraphic weaknesses. We envision that the metamorphic Veleta basement would behave more as a strong and coherent continental crustal block compared to unmetamorphosed sedimentary rocks. Additionally, even though some lithologies similar to the Veleta appear in the higher units, their intercalation with various lithologies such as metabasites and serpentinites result in weak and negatively buoyancy pockets that differ from the homogenous thick and continuous Veleta sequence.
Despite being subducted in the Cenozoic, our thermochronometers coupled with microstructures suggest that the Veleta basement did not exceed ∼450°C, which contradicts previous P/T estimates. Similar controversies are present in the Alps, where garnet and zircon ages are ∼320-340 Ma (Liati et al., 2009;Sandmann et al., 2014) and preserve Variscan metamorphism rather than later Alpine Orogeny metamorphism. This discrepancy merits a reevaluation of P/T estimates based on various methods or further investigation into the conditions under which minerals grow and reset their thermo-chronometer clocks during metamorphic events.

The Imbrication of the Sedimentary/Igneous Cover (Mulhacén Succession)
The Mulhacén succession forms a Permian to Mesozoic syn-and postrift sedimentary/igneous package with heterogeneous lithologies, including mostly metasedimentary rocks in addition to orthogneisses, metabasic rocks, marble and ultramafic slivers (Figure 10a; e.g., Poulaki & Stockli, 2022;Puga et al., 2005;Menzel et al., 2019). During Cenozoic subduction, the Mulhacén succession was cut by a series of imbricate thrust sheets and experienced significantly higher temperatures during HP/LT metamorphism than the underlying Veleta basement unit. In contrast to the monotonous Veleta basement, the heterogeneous architecture with varying lithologies in the Mulhacén succession may have contained presubduction inherited faults or a mechanical stratigraphy that allowed for the localization of strain and formation of imbricate thrusts to stack tectonic slivers decoupled from the intact basement. Additionally, as has been shown in other subduction zones, underplating of sediments is largely facilitated by the involvement of lithological heterogeneities (Tewksbury-Christle et al., 2021). Even though some of the lithologies in the Mulhacén succession are similar to Veleta, the rheological layering likely assists the development of large thrust faults, which also lead to different P/T conditions. This scenario would be similar to basement-detached imbrication during fold-thrust belt development (Pfiffner, 2016). Since Veleta and Mulhacén indeed have somewhat similar lithologies, we conclude that the primary factors controlling the different structural styles are the unit thickness, metamorphic history, and presubduction architecture.
These findings have important implications for dynamics within subduction systems since the basement plays an important role in the forces that drive and resist subduction, as well as the exhumation processes. As previously proposed for the Western Mediterranean, the main driver of exhumation for the Betic subduction complex is slab rollback (e.g., Brun & Faccenna, 2008;Gautier et al., 1999;Jolivet et al., 2003) with the contribution of buoyant continental crust and the basement rocks that today represent graphitic mica schist at cores of elongated extensional domes. The involvement of thick and buoyant basement rocks in the Mediterranean is one of the major drivers of underplating (Agard et al., 2009(Agard et al., , 2018Jolivet et al., 2005) and contributes to the mass transfer and recycling of older continental crust (Doin & Henry, 2001).

Subduction Zone Coherent Underplating
Underplating along the deep subduction interface has been imaged by geophysical studies along convergent margins such as Chile, North-West America, Alaska, Japan, and New Zealand, as summarized by Scholl (2021). All these studies imaged underplating at depths of ∼20-30 km, but due to limitations of the methods, crucial information is missing from the deeper underplating zone down to ∼60 km depth where we know these rocks experienced metamorphism before returning back to the surface. At shallower depths, underplating may be initiated as the decollement steps down and accretes slivers to the base of the overriding plate (e.g., Bangs et al., 2020). Furthermore, the involvement of elevated topography may provide physical mechanisms or enhanced vertical buoyancy forces to detach crust and sediments from the sinking slab (e.g., Cloos, 1993). Work by Piana Agostinetti and Faccenna (2018) shows imbrication and underplating occurring at the Adriatic continental convergent zone at about 40 km depth using active seismic data and receiver function analysis. In this study, we were able to identify contacts that represent fossil thrust imbricate faults that accommodate underplating and shed light on the nature of these contacts (Figures 10b and 10c). Seismic reflection data have also shown that the plate boundary often changes from a single bright reflection at shallow depths to a wide package of bright reflectivity, which may indicate underplating at greater depths (e.g., J. Li et al., 2015).
Our systematic high-density sampling and petrographic observations from the NFC show continuous stratigraphic successions and multiple older-on-younger relationships at the structurally higher levels. Although the formation of metamorphic zircon rims depends on many factors such as mineral budget and fluid composition, we can assess general trends in their presence or absence from our structural column to infer the metamorphic conditions and timing of imbrication. The absence of Eocene zircon rims in several samples from the hanging wall and different apatite ages across the thrust faults suggests that stacking and underplating began during or after peak metamorphism and continues while these rocks are in the exhumation path (Figure 10b). The timing of their formation is also supported by the recent study by Porkoláb et al. (2022), where in situ Ar/Ar analyses suggest that imbrication occurs after peak metamorphism with contraction continuing until the units pass through the brittle-ductile transition (Behr and Platt, 2012;Booth-Rea et al., 2005. This is further supported by the large distribution of temperatures that have been revealed from various studies for the NFC as recently compiled by Santamaría-López et al. (2019). The overall higher temperatures in the upper NFC have also been previously observed for these units (Behr & Platt, 2012;Santamaría-López et al., 2019). These studies proposed an inverted metamorphic gradient; however, the lower pressures in Veleta than in the Mulhacén succession Menzel et al., 2019), as well as the different metamorphic conditions between these units, instead argue that these observations can be reconciled by different maximum depths of subduction ( Figure 10). In this scenario, the Mulhacén succession would be subducted first and reach peak metamorphic conditions starting in the Eocene and prolonged until the Early Miocene. The heterogeneous layering of lithologies within the Mulhacén succession allowed for the localization of strain and formation of imbricate thrust sheets around peak metamorphic conditions near the decoupling depth. These thrust sheets continued to develop during exhumation as they were subsequently stacked onto the trailing Veleta unit at shallower depths in the Late Miocene. Together, the culmination of a thick, buoyant package drove underplating of the fused Veleta and Mulhacén succession to the overriding plate (Figure 10b). These processes may have varied slightly between ranges in the Betic Cordillera due to different depths of NFC decoupling indicated by the partially recrystallized apatite in Sierra de los Filabres and Sierra Alhamilla suggesting that they experienced lower temperatures than Sierra Nevada (Figures 5e and 5f). While underplating is often thought to involve the formation of mélange zones and chaotic mixtures (Cloos & Shreve, 1988;Shreve & Cloos, 1986), we present concrete evidence for coherent underplating of imbricate stacks of tectonic slivers thrusted on top of each other. In many cases, the base of these thrusts is characterized by orthogneiss slivers. The differences in material rheologies might be due to the weak decoupling horizons where fault formation occurs. Similar coherent structures have been observed in the Cyclades (Kotowski et al., 2022;Poulaki et al., 2019;Seman, 2016) and on Syros with mixed brittle-ductile deformation where brittle blueschist veins are coeval with ductile shearing (Kotowski & Behr, 2019;Kotowski et al., 2022). Underplating of the NFC likely occurred around peak metamorphic conditions at ∼40-60 km depth, approximately where the slab reaches the Moho depth of the overlying Iberian plate. Duplex structures like the ones we documented here have been imaged in active subduction zones (Calvert et al., 2011;Henrys et al., 2013;Moore et al., 1991) but are also predicted to occur during the underplating process by geodynamic models (e.g., Menant et al., 2019). In contrast to this coherent manner of underplating, subduction zones dominated by pelagic sediments atop oceanic crust can potentially form a chaotic mélange during underplating (e.g., Angiboust et al., 2013;Cloos & Shreve, 1988). However, recent work by Harvey et al. (2021) shows decoupling and coherent underplating beneath a mélange zone by juxtaposition with the amphibolite domain of the slab, suggesting that even in mélange-like underplating, mixing may occur in a more systematic manner than previously thought. Many studies have shown that slabs often stall and stagnate at the mantle transition zone as a possible mechanism for rollback and trench retreat (e.g., Agrusta et al., 2017;Billen, 2010;Marquardt & Miyagi, 2015;Torii & Yoshioka, 2007). Therefore, if the Betic slab could not easily penetrate into the lower mantle, this could be a mechanism for triggering slab rollback during Miocene, underplating of the NFC, and the massive trench migration toward the West as has been documented by various studies (e.g., Brun & Faccenna, 2008;Faccenna et al., 2004).
The presence of Miocene apatite and zircon ages with metasomatic morphologies located only within the upper NFC records metamorphism under the presence of abundant fluids. Hence, the previously formed thrust faults are not only potential zones of localized strain but also act as weak zones during the exhumation processes and pathways for fluids during various fluid pulses that have been identified in the upper units of the NFC (Figure 10c; Dyja et al., 2016). Additionally, Late Miocene extensional shear zones along the NFC units have been previously proposed by various authors, and the thrust faults identified in this study would be ideal inherited weak zones for subsequent reactivation during unroofing , 2010Porkoláb et al., 2022). Rapid rollback and trench migration could potentially facilitate slab tearing and give rise to the reheating, abundant fluid percolation, and punctuated arc magma genesis observed in the Late Miocene. Similar polyphase metamorphic relationships have been observed in other tectonic settings in the western United States (Constenius et al., 2003) as well as in Cyclades, Greece (Poulaki et al., 2019(Poulaki et al., , 2021.

Conclusions
In this study, zircon and apatite geo/thermochronology coupled with microstructural and field observations provide unprecedented details on the underplated subduction complex in Southern Spain. Within the NFC, we identify at least two phases of subduction zone metamorphism in the Eocene and Miocene as well as an increase in the metamorphic grade from the lower to upper units. Our findings indicate that subduction must have initiated at least since ∼60 Ma, which is significantly earlier than previous models. The polyphase metamorphic history of the NFC includes diachronous HP/LT peak metamorphic conditions in the Eocene and Early Miocene, followed by reheating and fluid-driven metasomatism during underplating and exhumation in the Late Miocene.
This work highlights the importance of the presubduction architecture for structural styles of deformation during underthrusting and underplating. The Veleta basement of the NFC had already experienced metamorphism during the Variscan orogeny and Jurassic rifting prior to subduction. Upon entering the subduction system, it behaved as a thick, buoyant block that resisted deep subduction and acted as an internal force mechanism to assist underplating and transfer to the overriding plate. In contrast, the structurally higher Mulhacén succession of the NFC formed coherent slivers with imbricate thrusts coeval with peak metamorphism above 500°C. Sedimentary protoliths and stratigraphic layering play an important role in where deformation becomes localized within thicker shear zones (∼100s m) and strain localization of thrust faults, which appear to preferentially develop along the boundary of more buoyant lithologies, such as slivers of orthogneisses. Decoupling of subducted materials from the downgoing slab occurs during or after peak metamorphism. For the Betic subduction zone, density contrasts between metamorphosed sediments and the overriding mantle wedge, the trailing subduction of a buoyant coherent block, and slab stagnation near the mantle transition zone may have been contributing factors to underplating the NFC, transitioning to exhumation, and trench migration to its current position at the Gibraltar Arc.

Data Availability Statement
All of the data are publicly available at https://doi.org/10.5281/zenodo.7644485. All the geochronology data can also be found at www.geochron.org.