Spatio-temporal patterns of Pyrenean exhumation revealed by inverse thermo-kinematic modeling of a large thermochronologic data set

The feedback between tectonically driven rock trajectories and erosion controls orogenic topography, generates exhumation, and integrates crustal deformation, lithology, surface processes, mantle dynamics, and climate (e.g., Whipple, 2009; Jamieson and Beaumont, 2013). Thermochronology is the prime tool for quantifying exhumation and, when paired with independent constraints on topography, can provide valuable information for investigating these relationships. A large spatial- and temporal-scale approach facilitates investigation into the rates, patterns, and drivers of orogen-wide exhumation within a regionally consistent framework. The Pyrenees mountains (France/Spain) present a unique study area because they comprise a small orogen with a well-constrained tectono-sedimentary evolution, for which spatially extensive thermochronologic data are available. Additionally, a recently published model of orogen-scale topographic evolution (Curry et al., 2019) constrains surface uplift through time, providing a rare opportunity to investigate spatio-temporal variations in both exhumation and rock uplift.

Pyrenean exhumation is variably attributed to (1) southward-propagating deformation of an orogenic wedge (Fitzgerald et al., 1999; Sinclair et al., 2005), (2) westward-propagating deformation due to variable pre-orogenic crustal extension (e.g., Puigdefàbregas et al., 1992; Whitchurch et al., 2011), and/or (3) base-level dynamics (Coney et al., 1996; Fillon and van der Beek, 2012). The localized nature of most Pyrenean thermochronologic studies prevents validation of a consistent model of exhumation; while several studies examine across-strike trends, few investigate along-strike variation. We use inverse thermo-kinematic modeling of a large, compiled multi-thermochronometer data set (n = 264) and incorporate topographic boundary constraints to evaluate controls on orogen-scale exhumation in the Pyrenees. We aim to explain two key features of the data: (1) an along-strike trend in exhumation, and (2) the existence of post-orogenic cooling ages. Our results show that spatio-temporal patterns of exhumation are controlled by structural inheritance and post-orogenic surface processes, and to a limited extent, by thrust-sheet kinematics.

The Pyrenean orogen formed due to Late Cretaceous to Miocene convergence between Europe and Iberia. Prior to convergence, both the Paleozoic Variscan orogeny and Mesozoic extension affected the region (e.g., Muñoz, 2019). Structural inheritance from this tectonic history persists in the modern orogen, with Pyrenean thrusts exploiting a now-inverted Mesozoic extensional fault system, producing the dominant east-west structural trend (Fig. 1). The mountain belt is divided into two oppositely verging thrust belts. The Southern Pyrenees thrust belt makes up the pro-wedge (downgoing plate) and is further divided into the Axial Zone (AZ) and the South Pyrenean zone (SPZ; Fig. 1). The AZ comprises an antiformal stack of Paleozoic basement thrust sheets in the eastern and central Pyrenees (Muñoz, 1992; Vergés et al., 2002). The SPZ is a south-vergent thin-skinned fold-thrust belt that propagates into the pro-foreland Ebro Basin (Fig. 1). The Ebro Basin became endorheic in the late Eocene, and >2 km of conglomeratic material was deposited and later excavated when the basin reopened in the Miocene (Coney et al., 1996). The North Pyrenean zone (NPZ) is a north-verging thick-skinned belt of inverted Mesozoic basins that makes up the retro-wedge (Muñoz, 2019).

We used the thermo-kinematic modeling code Pecube (Braun, 2003; Braun et al., 2012) to track rock-particle paths and temperatures during exhumation and to predict thermochronometer ages given an imposed tectonic (rock-uplift) and topographic (surface-uplift) scenario. We used an inversion approach to identify a tectonic scenario within a predefined parameter space that best fits the available thermochronometer data (details are provided in the Supplemental Material¹). Each inversion consists of 10,080 forward models that converge on a combination of input parameters resulting in the lowest calculated χ² misfit. We then assessed the inversion results by calculating the marginal posterior probability density function of each inverted parameter, thus providing a Bayesian estimate of parameter resolution (cf. Braun et al., 2012). Pecube modeling parameters are provided in Table S3 in the Supplemental Material.

Our goal was to determine whether the thermochronology data are best described by thrust sheet–controlled uplift (model A) or along-strike variable uplift (model B) (Fig. 1), where “uplift” implies rock uplift as defined by England and Molnar (1990). In model A, the AZ is divided into major thrust sheets that may move independently (Fig. 1C). In model B, the AZ is broken into four independent zones distributed east-west along strike (Fig. 1D). We treat the NPZ as distinct from the AZ due to its earlier exhumation history (Fitzgerald et al., 1999; Vacherat et al., 2016).

We invert for (1) rock-uplift rate through time, and (2) thickness of overburden deposited and subsequently removed in the Miocene (Fig. S1). In each inversion, the orogenic period (56–20 Ma) was divided into seven 5 m.y. time bins (earliest bin is 6 m.y. to coincide with the beginning of the Eocene at 56 Ma) to identify potential changes in rock-uplift rate (Fig. S1); for each time bin, the model inverts for the uplift rate that best reproduces the thermochronology data. The post-orogenic (<20 Ma) uplift rate was set at 0.01 km/m.y., following previous inversions (Fillon and van der Beek, 2012). The models include surface uplift (topographic evolution) based on the results of Curry et al. (2019), which involves topographic growth from 56 to 23 Ma and minor decay from 23 to 0 Ma (Fig. 1B). Including independently determined surface-uplift constraints is a novel approach that enables exploration of both rock uplift and exhumation.

We additionally ran two sets of inversions to test the influence of Miocene sediment blanketing on post-orogenic cooling patterns. In one set, we shut off rock uplift from 20 to 0 Ma but allowed as much as 2.5 km of deposition and erosion between 23 and 10 Ma (see the Supplemental Material for details). In another, we allowed tectonically driven rock uplift until 5 Ma but bar Miocene deposition.

We evaluated our models using 264 compiled thermochronologic ages, including apatite (U-Th)/He (AHe, n = 77), apatite fission-track (AFT, n = 123), zircon (U-Th)/He (ZHe, n = 40), and zircon fission-track (ZFT, n = 24) data (Fig. 2) (Morris et al., 1998; Fitzgerald et al., 1999; Sinclair et al., 2005; Gibson et al., 2007; Jolivet et al., 2007; Maurel et al., 2008; Gunnell et al., 2009; Metcalf et al., 2009; Herman et al., 2013; Bosch et al., 2016; Vacherat et al., 2016; Fillon et al., 2020). These data collectively record rock cooling from ∼300 to ∼60 °C and encompass the pre-, syn-, and post-orogenic periods. Ages vary from ca. 170 Ma for unreset ZFT in the NPZ to ca. 6 Ma for AHe in the western AZ; the bulk of the ages are between ca. 50 and ca. 20 Ma. The AFT and ZFT ages are single-sample, multi-grain “central ages”; the AHe and ZHe ages are either aliquot or average single-grain ages. All sample ages, uncertainties, locations, and references are given in Table S2.

Model B (along-strike variable uplift; Fig. 2) presents the best fit to the data, with a χ²-misfit value of 5.6; model A (thrust sheet–controlled uplift) has a misfit of 11.6 (Fig. S2). Inversion of the NPZ data produces a misfit of 4.9. The best-fit results of model B illustrate a westward-propagating wave of exhumation in the AZ during the orogenic period (Fig. 2). Until 40 Ma, the eastern AZ (zones 3 and 4) underwent moderate rock uplift, peaking at >1 km/m.y. between 35 and 30 Ma, after which uplift slowed to near zero (Fig. 2E). In the western AZ (zones 1 and 2), rapid rock uplift did not begin until 30 Ma, before which rates did not exceed 0.4 km/m.y. (Fig. 2E). Results for the NPZ show peak rock-uplift rates (0.7 ± 0.06 km/m.y.) between 50 and 40 Ma; otherwise, rates are <0.4 km/m.y. (Fig. S2). Uncertainties on the topographic evolution add 0.007–0.04 km/m.y. of uncertainty to these results. Although the thrust sheet–controlled uplift model does not fit the data best, examination of the patterns of cooling and misfit reveal a southward-younging trend in the east-central region of the AZ (zone 3; Fig. 3A).

The “sediment-blanketing” and “tectonic” models for post-orogenic cooling produce similar values of misfit (5.6 and 6, respectively), but the “tectonic” scenario entails a surge of rock uplift (∼0.5 km/m.y.) since 10 Ma in the western AZ. In the “sediment-blanketing” scenario, the thickness of deposited and subsequently eroded sediment favored by the best-fit model varies along strike, from >2 km in the western AZ to 0.4–1.2 km in the eastern AZ (Fig. 2E).

Using rock-uplift rates from the best-fit model and the topographic evolution, we calculated erosion of the Pyrenees through time (Fig. 3B). Results show that the eastern AZ underwent the most syn-orogenic erosion, but the western AZ underwent the most post-orogenic erosion (Fig. 3C). Prior to 30 Ma, the eastern AZ was eroded by 13–17 km and the western AZ by only 5–6 km. Between 30 and 20 Ma, the east eroded by 0.6–4 km and the west by ∼10 km; following 20 Ma, the western AZ was eroded by ∼2.6 km compared to 0.5–1 km in the eastern AZ (Fig. 3C). The final, cumulative magnitudes of erosion are nearly equivalent along strike at ∼17–19 km.

Based on our inversions, we infer a dominant first-order trend of westward-propagating rock uplift for the Pyrenees (Fig. 2E). We propose that the pre-orogenic tectonic template is responsible for the westward-younging rock uplift and exhumation pattern. Prior to convergence, the western Pyrenees comprised a wide extensional zone with a broader region of transitional crust than the eastern Pyrenees (Jammes et al., 2014; Teixell et al., 2016; Muñoz, 2019; Fig. S4). Hence, although the onset of shortening was synchronous along strike (Teixell et al., 2018; Muñoz, 2019), thinned crust was consumed earlier in the east than the west due to this inherited crustal configuration, leading to diachronous crustal thickening and uplift (Figs. 4A and 4B). Previous work proposed westward propagation of uplift and exhumation based on sedimentation patterns in the southern foreland basin, which show coeval marine deposition in the west and continental deposition in the east throughout the Cenozoic (Puigdefàbregas et al., 1992; Garcés et al., 2019). The southern fold-thrust belt also shows westward-younging deformation: the eastern frontal structures were active until the mid-Oligocene, whereas the western Sierras Exteriores were active into the Miocene (Fig. 4; Burbank et al., 1992; Millán Garrido et al., 1995). This thin-skinned frontal deformation is spatially and temporally linked to the rapid hinterland exhumation inferred here.

Several thermochronology studies around the ECORS (Etude Continentale et Océanique par Réflexion et Réfraction Sismique) seismic cross-section in the central Pyrenees inferred southward propagation of exhumation (Fitzgerald et al., 1999; Sinclair et al., 2005; Metcalf et al., 2009). We observe a similar trend in the area of those studies (Fig. 3A) but find no evidence for such a pattern along strike in the AZ, suggesting the findings from this region cannot be extrapolated to the rest of the Pyrenees. We speculate that this pattern is related to the stronger development of the antiformal stack in the central Pyrenees (our zone 3) compared to the western Pyrenees (Fig. S4; Vergés et al., 2002; Sinclair et al., 2005; Teixell et al., 2016).

An enigmatic aspect of the thermochronology data is the presence of middle to late Miocene (15–6 Ma; Fig. 2D) AHe ages in the western AZ. These post-orogenic ages are variably interpreted as recording renewed deformation (Jolivet et al., 2007; Bosch et al., 2016), valley excavation following the connection of the foreland Ebro Basin to the Mediterranean (Coney et al., 1996; Fitzgerald et al., 1999; Fillon and van der Beek, 2012), or enhanced erosion resulting from Pliocene climate change, as proposed for the European Alps (Cederbom et al., 2004). All interpretations are problematic. Contrary to the first interpretation, there is no structural or stratigraphic evidence for ongoing, late Miocene deformation in the AZ (Muñoz, 2019). The presence of young cooling ages both north of the drainage divide and in regions where the maximum thickness of the southern conglomerate is interpreted to be <2 km disputes the second interpretation. Finally, the along-strike variability of this cooling signal is inconsistent with a large-scale climatic influence.

We interpret these ages as recording post-orogenic erosion based on the lack of evidence for ongoing deformation and the ability of our model to replicate these cooling ages without tectonically driven uplift (Fig. 2D). We propose that three synergistic circumstances drove post-orogenic erosion: greater stream power due to base-level fall, isostatic rebound, and variable substrate erodibility. The Ebro Basin was endorheic from ca. 35 to ca. 10 Ma (Fig. 4; Garcia-Castellanos et al., 2003; Fillon and van der Beek, 2012; Garcés et al., 2019), during which time the basin was buried and the southern Pyrenees were backfilled with 1–3 km of conglomeratic material shed from the Pyrenees (Fig. 4C; Coney et al., 1996). When the Ebro river connected with the Mediterranean in the late Miocene, the ensuing base-level drop initiated a wave of erosion that propagated through the basin and removed hundreds of meters of sediment, prompting as much as 600 m of isostatic rebound (Garcia-Castellanos and Larrasoaña, 2015).

Bernard et al. (2019) argued that lithological variations impose a strong control on post-orogenic erosion in the Pyrenees. Mapping of the AZ delineates a boundary between rocks of the metamorphosed Variscan hinterland in the east and unmetamorphosed foreland rocks in the west (Fig. 4; García-Sansegundo et al., 2011). When the erosive wave reached the eastern AZ, it encountered more durable Variscan hinterland rocks that had undergone >18 km of exhumation during the Pyrenean orogeny (Figs. 3B and 3C). In contrast, in the western AZ, the erosive wave encountered more erodible rocks in a region that had undergone less Pyrenean exhumation (Fig. 3C), facilitating the 2–3 km of erosion necessary to exhume rocks from below the AHe closure temperature (Fig. 4).

Eroding 2–3 km of rock from the western AZ would have resulted in a ∼500 m reduction in topography, considering flexural isostasy and depending on crustal density. Taking into account the additional rebound associated with excavation of the Ebro Basin (∼100 m in the western Pyrenees; Garcia-Castellanos and Larrasoaña, 2015), this estimate is within the range of topographic decay (300–600 m) since the Oligocene proposed by Curry et al. (2019).

These results highlight the value of modeling large thermochronology data sets and the novelty of incorporating paleotopographic constraints. We propose a regionally consistent model to explain exhumation patterns in the Pyrenees, finding that the data are explained by a combination of pre-orogenic tectonic inheritance, syn-orogenic wedge dynamics, and post-orogenic surface processes.

This research is funded by TOTAL through the Coupling Lithospheric Deformation and Surface Processes (COLORS) project. H. Sinclair, D. Barbeau, and an anonymous reviewer provided helpful reviews.