The present-day topography in Iberia is related to geodynamic processes dealing with lithospheric-scale deformation. However, little attention has been paid to the role of inherited crustal- or lithospheric-scale structures involved in the recent observed large-scale topographic patterns. Whereas the analysis of brittle structures focuses on the evolution of Mesozoic sedimentary basins and their subsequent response to tectonic inversion, their contribution to mountain building has been underestimated. Large numbers of structures, from ductile to brittle, which affected the whole lithosphere, were developed during the evolution of the Cantabrian orocline (ca. 310–300 Ma). The contribution of these Paleozoic post-Variscan structures, together with lithospheric mantle evolution and replacement during orocline development in the Mesozoic and Cenozoic geological evolution of Iberia, remains unexplored. To explore the role of these inherited structures on the final configuration of topography during N-S Pyrenean shortening, we carried out a series of analogue experiments complemented by surface velocity field analyses. Our experiments indicate that strain was concentrated along preexisting crustal- to lithospheric-scale discontinuities, and they show several reactivation events marked by differences in the velocity vector field. Differences in fault displacement were also observed in the models depending upon preexisting fault trends. The obtained results may explain the different amount of displacement observed during the reactivation of some of the post-orocline structures in Iberia during the Cenozoic, indicating the key role of unveiled structures, which probably have accommodated most of the Alpine shortening.

Since the advent of the Wilson cycle concept (Wilson, 1968), reactivation of previous structures represents one of the main controls in the tectonic evolution of continents. When initiating a Wilson cycle by opening new oceans, the reactivation of previous suture zones is likely to nucleate the initiation of new oceanic realms (e.g., Burke et al., 1976; Bailey et al., 2000; Tikoff et al. 2001; Murphy et al., 2006; etc.). Rift initiation through fault reactivation is better understood in segments of orogens depicting a rather linear attitude (Thomas, 2006; Murphy et al., 2006). However, in complex curved orogens, as in the case of the western European Variscan belt, it is still difficult to ascertain (Fig. 1A).

The tectonic evolution of north-central Iberia since the latest Carboniferous includes several deformation episodes postdating the Variscan orogeny: (1) formation of a curved orogen, known as Cantabrian orocline (i.e., Gutiérrez-Alonso et al., 2004, 2008, 2012, 2015; Weil et al., 2010, 2013; Martínez-Catalán, 2012; Pastor-Galán et al., 2016, 2017; Fernández-Lozano et al., 2016; Murphy et al., 2016); (2) opening of the Mesozoic rift-related Bay of Biscay (Sibuet and Collette, 1991; García-Mondéjar, 1996); and (3) the Alpine convergence history (among others, Cloetingh et al., 2002; Vergés and Fernández, 2006; Casas-Sainz and de Vicente, 2009; de Vicente and Vegas, 2009).

However, not much attention has been paid to the putative changes produced in the Iberian lithospheric mantle during orocline development (i.e., Gutiérrez-Alonso et al., 2011a, 2011b), and the meaning and origin of late orocline structures attributed to the tightening of the 180° bend that defines the Cantabrian orocline. These structures, including faults observable at the present-day erosion level, controlled the basement structural grain in northern Iberia (Fig. 1A), and their imprint is noteworthy on the subsequent Mesozoic and Cenozoic belts and basins development, on the Alpine evolution of the northern and central Iberia mountain chains, and probably on the present-day relief of western continental Europe.

One of the consequences of generating 180° curved lithospheric-scale oroclinal bends is that they are not able to be bent any further, and if the strain conditions that caused the orocline persist, additional shortening is assumed by tectonic structures. Previous analogue experiments carried out by Pastor-Galán et al. (2012) suggested that conjugated strike-slip faults (shear zones in depth), which crosscut both parallel limbs of the orocline, develop shortening normal to the orocline axial plane and extension parallel to it (Fig. 1A). Similar lithospheric-scale fault patterns have been described related to the formation of the Himalayas, where lithospheric shortening produced conjugate faults that affected the Eurasian plate. These structures accommodate part of the general N-S–trending shortening caused by the collision of Eurasia with the Indian continent and dividing the Tibet, Tarim, and Tien-Shan realms (i.e., Tapponnier and Molnar, 1979; Tapponnier et al., 1987; Avouac and Tapponnier, 1993; Thatcher, 2007; England and Houseman, 1986; England and Molnar, 2005).

In Iberia, post-Variscan structures have long been described to play a key role during Mesozoic and Cenozoic time. However, difficulties in the recognition of the Alpine evolution and role of the preexisting faults have led several authors to suggest different mechanisms for intraplate mountain building, including: (1) the westward transmission of stress along NW-SE–trending structures (Martín-González, 2009; Martin-González and Heredia, 2011); (2) effective transmission of stress into the plate interior along NE-SW and NW-SE large-scale strike-slip faults (de Vicente and Vegas, 2009); and (3) strain partitioning mechanisms responsible for the observed complex pattern of topography (de Vicente et al., 2009, 2018).

This paper explores the influence of the alleged post-orocline structures together with the thermomechanical variations of the mantle lithosphere caused by the Cantabrian orocline development (Gutiérrez-Alonso et al., 2011a, 2011b). We investigated the contribution of preexisting structures to the present-day topography and the tectonic influence of strain transmission during the N-S Pyrenean orogeny, caused by the relative movement between Africa and Iberia. In order to understand the effects of the possible lateral strength variations caused by the putative differences in the post–Cantabrian orocline lithosphere underlying Iberia, we propose a kinematic model based on an integrated approach that combines analogue modeling, particle image velocimetry (PIV) analysis, and regional geology. Our results shed light on the evolution and kinematics of late orocline structures and lithospheric mantle variations and their contribution to recent deformations and present-day topography in Iberia.

During the final tightening of the Cantabrian orocline, in Permian times, the Variscan lithosphere was fractured into a conjugate system of NE-SW– and NNW-SSE– to NW-SE–trending strike-slip faults (Fig. 1A). This fracture pattern can be identified across central and eastern Iberia and in the French Central Massif (Fig. 1A; Vegas, 1975; Faure, 1995). In Iberia, this conjugate system partially controlled Permian, Triassic, and Mesozoic sedimentation and, subsequently, in many cases, like in the Iberian Chain, localized the position of the Alpine deformation fronts (Arthaud, 1975; Vegas, 1975; Casas-Sáinz and Gil-Imaz, 1994; Doblas et al., 1994; Guimerà et al., 2004; de Vicente et al., 2009). These crustal-scale structures are interpreted to have had horizontal displacements of ∼35–50 km, and, in some cases (Plasencia fault, Vilariça-Bragança fault system, Lower Tagus fault zone), they have reached deep into the lithospheric mantle, leading to the injection of Triassic-Jurassic boundary–age Central Atlantic magmatic province mantle-derived mafic melts (Ribeiro et al., 1990; Vegas, 2000).

According to the different characteristics of the Mesozoic geological evolution of Iberia, two main domains can be defined: The western Iberian Massif, which is devoid of large Mesozoic basins (except for the Lusitanian Basin, the origin of which is related to the opening of the Atlantic Ocean), and the eastern Mesozoic basins (Fig. 1B). The large conjugated strike-slip faults generated after Cantabrian orocline development were reactivated during the Mesozoic and Cenozoic with different kinematic regimes in both regions (Arthaud, 1975; Vegas, 1975; de Vicente et al., 2009). The different geodynamic evolution of these areas during this period has contributed to large lithospheric rheological variations. These thermomechanical differences represent a boundary between a somewhat stable Mesozoic western Iberian Massif, characterized by relatively cold and stable lithosphere, and a thinner, warmer, and younger lithosphere along the extended eastern sector of Iberia (Fig. 1B).

It is widely accepted that these thermomechanical differences between the western Iberian Massif and the eastern Mesozoic basins may have influenced the reactivation and tectonic inversion of late orocline–related faults (de Vicente and Vegas, 2009; Fernández-Lozano et al., 2012). Both domains are described in the following sections, together with a summary of their Alpine (Cenozoic) evolution.

Alpine Tectonic Evolution

The early episodes of the Alpine (Pyrenean) orogeny caused the deformation of the Iberian Massif and inversion of the NE basins and subsequent uplift during the Late Cretaceous–Eocene and the Oligocene–early Miocene, respectively (Casas-Sainz, 1993; Gómez et al., 2002; Guimerà et al., 2004; de Vicente et al., 2009; Teixell et al., 2018). As a result, the main Duero and Ebro continental foreland basins and the Tagus intraplate basin emerged, shaping the current geography of Iberia. Because of the different evolution of the Iberian Massif and the eastern Mesozoic domain, the resulting structures in both of them are significantly different.

Iberian Massif

N-S shortening was responsible for the building of the Cantabrian Mountains and the uplift of Paleozoic basement along the Iberian Massif and in the Central System (Spanish and Portuguese Central System), with shortening of ∼90 km (Alonso et al., 1996; Gallastegui et al., 2002; de Vicente and Vegas, 2009; Tavani and Granado, 2015; Llana-Fúnez and López-Fernández, 2015; Quintana et al., 2015; de Vicente el al., 2018), and the Cantabrian continental margin, decreasing progressively to the west (Gallastegui, 2000; Gallastegui et al., 2002; Pedreira et al., 2007; Fernández et al., 2015) into predominantly strike-slip displacements (Santanach et al., 1988; Santanach, 1994). These differences are highlighted by strong changes in the depth to the Moho, from 28–30 km offshore to 34–41 km beneath the main mountain uplifts (Montes de León, Spanish Central System, and Cantabrian Mountains, respectively; Pedreira et al., 2007; Díaz and Gallart, 2009; Díaz et al., 2015; Torne et al., 2015). Deformation started in the late Eocene and was well recorded in the Cantabrian continental shelf (Pedreira et al., 2015; Cadenas et al., 2018, and references therein). The cause of this shortening was the underplating of part of the Bay of Biscay crust under Iberian crust, as revealed by seismic data (Alvarez-Marrón et al., 1996, 1997; Gallart et al., 1995; Fernández-Viejo et al., 1998, 2000, 2012; Ayarza et al., 2004; Gallastegui, 2000; Gallastegui et al., 2002; Pedreira et al., 2003, 2007; Pulgar et al., 1996).

Strike-slip fault activity has been broadly constrained for the Cenozoic (mainly Eocene–Oligocene–early Miocene activity) through the use of fission tracks (see above) and ages of syntectonic deposits (de Bruijne and Andriessen, 2000, 2002; Barbero et al., 2005; Martín-González, 2006; Martín-González et al., 2008, 2012; Fillon et al., 2012; de Vicente and Muñoz-Martín, 2012; Fillon et al., 2016). In addition, fault reactivation in western Iberia is suggested by the presence of structural and geomorphological features, including paleostress indicators and the rapid fluvial incision observed since the Miocene in some areas (Antón et al., 2010, 2012). The NNW-SSE fault systems are mainly located in the northwestern part of Galicia, and they are characterized by right-lateral displacements concentrating present-day seismic activity in the western corner of Iberia (Santanach, 1994; Llana-Fúnez and López-Fernández, 2015). This system includes structures such as the As Pontes fault system, which shows compressional step-overs leading to the formation of depressions (Ferrús i Pinyol et al., 2005).

Eastern Mesozoic Basins

The observed differences in tectonic inversion style and the presence of a thick Mesozoic sedimentary cover deposited during an episode of rifting (>8 km of Jurassic and Cretaceous sediments accumulated along E-W basin depocenters; Hernando-Costa, 1973; Alvaro et al., 1979; Raven and van der Pluijm, 1986; Sopeña et al., 1988; Gómez-Pérez et al., 1998; López-Gómez et al., 2002; Aurell et al., 2002; Martín-Chivelet et al., 2002; García-Mondéjar et al., 2005; Omodeo-Salé et al., 2014, 2017) were controlled by the thermomechanical properties of the lithosphere. Crustal extension favored thermal subsidence and established the final configuration of crustal blocks limited by faults across the Basque-Cantabrian and the Iberian Basins.

Present-day crustal thickness ranges ∼36–44 km across the inverted Basque-Cantabrian and Iberian Basins, with a maximum up to 50 km in the central-west Cantabrian area, indicating a strong correlation between uplifted areas and Moho depth (Carballo et al., 2015; Díaz et al., 2015; Guimerà et al., 2016; Torne et al., 2015; Mancilla and Diaz, 2015). The inversion history of the basins started in the late Eocene (Guimerà and Alvaro, 1990; Solo et al., 2011), driven by the large-scale geodynamic configuration established during the Cenozoic convergence between Europe and Africa. N-S compression led to uplift and crustal thickening, and the formation of surrounding foreland basins such as the Duero and Ebro Basins (Millán et al., 1995; Casas-Sáinz and Maestro-González, 1996; Alonso et al., 1996; Casas-Sáinz et al., 2000; Gaspar-Escribano et al., 2001; Suárez-González, 2015). Evidence of Oligocene–early Miocene uplift was recorded by the formation of planation surfaces, which were subsequently captured during the establishment of the present-day fluvial network (Casas-Sáinz and Cortés-Gracia, 2002).

Cenozoic tectonic inversion of the Iberian Basin was associated with strike-slip displacement along NW-SE structures with compressional step-overs, which may have configured the strain partitioning scenario responsible for the mountain uplifts reported in the Iberian Range (de Vicente et al., 2009). In addition, observed E-W thrusts are responsible for fold wavelengths up to 5–13 km affecting the Variscan basement (Guimerà et al., 2004), and NW-SE– and NE-SW–trending faults maintain present-day shallow seismicity (<15 km deep), indicating neotectonic activity.

To explore the role of preexisting Variscan structures during the N-S Alpine shortening, we performed a study based on analogue experiments and surface analyses (i.e., analysis of particle displacement). Surface analyses consisted of a series of images acquired during modeling that were processed to obtain a surface velocimetry model from particle image velocimetry (PIV) analysis (Leever et al., 2011). The results re-create the fault evolution and distribution of strain during the main Alpine intraplate mountain-building episode that gave rise to the present-day topography in Iberia.

Analogue Modeling

These experiments complement previous work carried out by Fernández-Lozano et al. (2011) and Pastor-Galán et al. (2012), wherein they studied the role of shortening during the Alpine N-S shortening phase and after the closure of the Cantabrian orocline. The role of tectonic structures was tested by comparison of two sets of analogue experiments, where strips of weak silicone simulated the preexisting post-Variscan orocline–related fault zones.

The rheological properties of the lithosphere models were based on existing geological and geophysical data from Iberia. Models consisted of two different setups, where rheological differences were tested according to constraints provided by previous field studies (Fig. 2A; Suriñach and Vegas, 1988; Guimerà et al., 1996; Tejero and Ruiz, 2002; de Vicente and Vegas, 2009; Díaz and Gallart, 2009; Jimenez-Diaz et al., 2012; Carballo et al., 2015; Seillé et al., 2015; Torne et al., 2015). The experiments consisted of three layers, characterized by: a brittle upper crust, a ductile lower crust, and an upper lithospheric mantle (Fig. 2B). These layers rested over a high-density asthenospheric fluid inside a Plexiglas tank. A single moving wall deformed the model by 20% of bulk shortening according to rates of shortening suggested by de Vicente et al. (1996), and the side walls were lubricated to avoid side effects. Modeling materials were designed to represent two different types of lithosphere, a relatively hotter and younger lithosphere to the east (Iberian Basin), and a colder and older lithosphere to the west, highlighting the rheological differences in the lithosphere prior to Alpine shortening, i.e., differences that could have been inherited from post-Variscan orocline–triggered lithospheric thinning and delamination and intensely modified in the east due to Mesozoic extension. Details on material properties and scaling parameters are given in Table 1.

In addition to the simplest scenario (model A), characterized by a single weak zone representing a late Variscan structure, we compared the results with a more complex setup implemented with several weak zones (purple silicone strips in Fig. 2) representing the main E-W Late Cretaceous depocenters and large-scale shear zones (model B; de Vicente and Vegas, 2009). Similar setups have been previously implemented to study the effect of inherited heterogeneities in lithosphere-scale analogue experiments by Calignano et al. (2015a, 2015b, 2017). The role of the E-W–trending batholith in central Spain was also investigated according to the effect caused by topography and the rheological differences between the metamorphic basement and the Central System Granites (Martín-Velázquez and de Vicente, 2012). This rheological variation was modeled using a strip of weak silicone to simulate heat transfer from the granite body (see Table 1).

The physical properties and scaling parameters of the models followed previous studies carried out by Fernández-Lozano et al. (2011, 2012) and Sokoutis and Willingshofer (2011), in agreement with geometric and dynamic similarities proposed by Ramberg (1967b) and Weijermars and Schmeling (1986). The geometric similarity was achieved by assuming the following relationship according to scaling properties suggested by Brun (1999):
where σ refers to stress, ρ is density, g is gravitational acceleration, and L is the length scale. The asterisk refers to the ratio between model and nature. Since g* = 1 in our study, and the densities of silicon putties (∼1400–1500 kg m−3), and a wide variety of rocks in nature (2600–3000 kg m−3) are in the same order of magnitude, we can assume that σ* = L*. Therefore, considering the physical material properties and their equivalents in nature, the stress ratio is 3.25 × 10−7, which implies a geometric scaling of 1 cm in the model to ∼15 km in nature.

In addition, dynamic similarity was obtained through the dimensional analysis proposed by Ramberg (1967a) and Weijermars and Schmeling (1986), based on the Ramberg (Rm) and Smoluchowsky (Sm) numbers (see Sokoutis et al., 2007).

PIV Analysis

The analysis of fault kinematics carried out on the model surface was implemented following the methodology described in Leever et al. (2011) and using the open source software MatPIV (Sveen, 2004). Their method is based on study of the particle displacement field performed over the surface of the experiments. The PIV method compares two images within a fixed time interval. The image evaluation is carried out by dividing the PIV recording into several small subareas called “interrogation windows.” The record is measured in picture elements or pixels. Therefore, if small-grain particles are scattered over the model surface, those particles within the interrogation windows will be correlated.

Evaluation of PIV recordings was conducted using the cross-correlation method, defined as a standard statistical method of estimating the degree to which two compared series of data are correlated (Westerweel, 1997; Raffel et al., 1998). The PIV method aimed to investigate fault slip by comparing the velocity field and slip orientation during model deformation. The method is based on the calculation of a directional derivative, δυ, from the original velocity field by subtracting adjacent vectors with different orientations,
where υ is the initial vector field, and δυ is its directional derivative. Consequently, if the displacement field is constant, the subtraction of adjacent vectors provides no resultant. However, if the magnitude and orientation of the vector field change, it will provide a different vector ≠ 0. When the length of one of the vectors in δυ exceeds a threshold value, its azimuth is plotted in a color scale (0–360°). Therefore, high angles indicate convergence directions (∼90°), while low angles (<45°) show dextral strike slip. Following this reasoning, we were able to identify different mechanisms of fault slip (oblique or parallel to convergence direction) that took place during model deformation.

In this way, this method allowed us to study and distinguish among active, reactivated, and nonactive faults. That means that we were able to study regions inside the models where strain localization was taking place during deformation. Moreover, the method provided new insights into the strain distribution and mechanism of partitioning inferred from the particle velocity field.

Two sets of lithospheric-scale model runs were performed under a unidirectional convergence direction. Material properties and scaling parameters are summarized in Table 1 (Fig. 2B).

Model A

The rapid evolution of deformation was controlled by the presence of preexisting weak zones and the net discontinuity between both lithospheric types, localizing the strain distribution. The weak zones were immediately reactivated, showing two stages of deformation related to strike-slip displacements followed by thrusting producing between 10% and 15% of bulk shortening (Fig. 3). The surface velocity field calculated with the PIV method implemented by Leever et al. (2011) revealed counterclockwise rotation of velocity vectors caused by displacement along the main fault limiting the two lithospheric types. However, clockwise rotation of slip vectors was also recorded inside the stiffest lithosphere, along and close to the preexisting central weak zone, which was at a large angle to the boundary between the two different types of lithospheric mantle. These changes observed in vector distribution support the notion of different reactivation of crustal heterogeneities during compression according to their initial orientation to the prevailing stress field.

The fault displacement regime was affected by subsequent episodes of reactivation, changing from pure strike-slip to thrust components and leading to topographic uplift (between 5% and 20% of bulk shortening). Back-thrusting occurred at 15% and 20% of bulk shortening in the central part of the model. In addition, non-coeval episodes of thrusting (fore- and back thrust) led to compensation of the observed uplift through a series of triangle zones that accommodated the amount of displacement caused by thrust reactivation. Strike-slip movement was apparently restricted to the main oblique discontinuities deformed under a regime of transpression (i.e., lithosphere boundaries and weak zones). Consequently, a component of transpression was always associated with fault reactivation (compare top-view images with digital elevation models in Fig. 3).

The modeling results show that faulting resulted mostly in imbricated thrusts and back thrusts producing pop-up structures, whereas wide, flat, and highly elevated regions formed between the main thrust systems. Distribution of brittle deformation was focused within strong lithosphere. However, deformation of the weak lithosphere resulted in strain localization along closely spaced thrust systems. It is worth noticing the interference deformation patterns among E-W, NE-SW, and NW-SE crustal structures occurring under unidirectional shortening (coexisting strike-slip and thrusting components), suggesting strain partitioning limited by the previous faults, either crustal or lithospheric. Although strain partitioning seems to have been strongly influenced by the presence of lateral strength variations affecting the lithosphere, the initial position of preexisting weak zones also appears to have played a major role in the mode through which strain was spread over the surface of the model. These weak zones also defined the boundary of different strain transmission domains.

Clear differences in the topographic elevation above the different types of lithosphere are highlighted by the digital elevation models (DEM) shown in Figure 3. The weak lithosphere shows high uplift rates and deep depressions (compare DEMs in Fig. 3). These differences are also observed in the relative magnitude of displacement calculated on the model’s surface. The central part of the model represents a major boundary for displacement vectors, probably associated with the main preexisting zones of weakness (5%, 10%, 15% bulk shortening). Major rates of uplift were found at 20% bulk shortening along the weak lithosphere, as indicated by the position of uplifted areas.

Model B

Model B was implemented with four weak zones representing the main Iberian margin conjugated fault system (NW-SE faults comprising the Plasencia, Vilariça-Bragança-Manteigas, and Penacova-Regua-Verin faults and the Lower Tagus fault zone; and the NE-SW system, including the Somolinos, Cantabrian, Ubierna, Pamplona, Demanda, and Altomira faults). Moreover, a weak zone within the model upper crust was used to represent the granitic belt that entirely crosses the Spanish Central System (oriented around 80° with the direction of convergence). It is also important to remark that the main E-W Cretaceous depocenters were also modeled using strips of weak silicone in the upper mantle.

The evolution of the topography in the model was strongly influenced by fault growth. Tectonic thrust structures caused the elevation of the model surface at successive steps (Fig. 4; compare DEMs). Rapid uplift of the model surface was caused first inside the weak lithosphere, following the orientation of the main (E-W) depocenters and the boundary between both types of lithosphere (trending NW-SE) at ∼10% bulk shortening. Uplift along the modeled depocenters was caused by thrusting, whereas dextral strike slip across the lithosphere boundary involved a component of transpression. The analysis of surface displacement of particles indicated a counterclockwise rotation of the vectors along the southern tip of the lithosphere boundary and the presence of oblique velocity vectors around the entire discontinuity. Moreover, the PIV analysis also recorded movement along the preexisting weak zones trending NW-SE, NE-SW, and E-W. No deformation was recorded along structures oriented NNE-SSW within the strong lithosphere. Successive episodes of displacement were observed along NE-SW structures after reactivation during further deformation, finally leading to left-lateral strike-slip movements combined with thrusting (uplift between 5% and 10% bulk shortening). Similarly, the weak zone representing the Spanish-Portuguese Central System granitoids caused an abrupt uplift along two well-oriented E-W thrusts facing north and south, and defining a crustal pop-up structure. This part was reactivated at 10%, 15%, and 20% of bulk shortening (see vectors in Fig. 4). The thrusts’ vergence significantly changed from one lithosphere to another. Generally, north-verging thrusts were characteristic of a weak lithosphere (i.e., facing the moving wall), whereas the strong lithosphere was associated mainly with south-vergent thrusts and pop-ups, as was revealed from displacement vectors (see 15% bulk shortening in Fig. 4). Magnitudes of displacement vectors showed small differences between deformed areas. However, at 15% bulk shortening, a major rate of displacement was observed to occur along the weak lithosphere, probably enhanced by reactivation of NW-SE weak zones. This situation suddenly changed, reaching higher-elevation magnitudes along the moving wall within the strong lithosphere sector, and it is interpreted to be a border effect.

The Mesozoic–Cenozoic tectonic evolution of Cantabrian orocline–related structures remains controversial and varies from east to west Iberia according to variations in lithospheric strength attributed to inheritance from the Mesozoic rifting episode (de Vicente and Vegas, 2009) together with changes in the mantle derived from orocline-triggered lithospheric delamination (Gutiérrez-Alonso et al., 2004). The tectonic inversion style in both regions also differs and was governed by the geometry and location of preexisting conjugate fault sets and the position of depocenters developed during the Mesozoic N-S extension that acted as buttresses (Fig. 5A). Altogether, we propose a new model that represents an analogue of other Alpine settings, including the Himalayas or the Alps (Peltzer and Tapponnier, 1988; Luth et al., 2013), where N-S shortening was accommodated in E-W reverse fault corridors, whereas the conjugate fault system represented by NE-SW to NNE-SSW and NW-SE structures played a scarp-like role as a result of strain partitioning between NE-SW strike-slip fault corridors and E-W thrusts (Fig. 5B). This model of distributed deformation in two sectors of Iberia with different lithospheric characteristics best fits the observed differences between the western Iberian Massif and the eastern Mesozoic basins. While other models based on a large-scale crustal detachment or a combination of lithospheric folding and strain partitioning have been proposed, none of them successfully and independently explain the differences in geometry and tectonic style observed between the Iberian Massif and the eastern Iberian Chain (Warburton and Alvarez, 1989; Cloetingh et al., 2002; Fernández-Lozano et al., 2012; Quintana et al., 2015).

However, our analogue experiments have shown that under N-S shortening, NE-SW structures showed sinistral strike-slip movements. These faults were successively reactivated, and in some cases developed a thrust component. The small displacement shown by NNE-SSW structures could be the result of a distribution of deformation toward the nearby Variscan structures, which was enhanced by the rotation into a NE-SW direction during the development of the Variscan orocline.

Despite the relatively large amount of displacement along NE-SW structures in the analogue experiments, geological evidence in central Spain has shown that displacement rates are low (i.e., <3–5 km for the 550 km Plasencia fault; Villamor-Pérez et al., 2012). These differences might be caused by the presence of minor thrusts affecting the Paleozoic basement next to the fault, which may absorb part of the deformation, or they may be the result of movements in opposite directions along the fault planes that have been obliterated. All in all, the vector field obtained over the model’s surface suggests a large deformation area with limited uplift along these structures when compared with the observations made along NW-SE faults (compare DEMs for models A and B in Figs. 3 and 4).

The heterogeneous crust in the central part of the experiments (weak zone representing granites from the Spanish and Portuguese Central System [SPCS] in model B; Fig. 4) was efficiently reactivated, leading to strain localization and subsequent tectonic uplift (see PIV results in Fig. 4 and DEM results). This weak zone controlled the evolution of the Spanish and Portuguese Central System range during deformation. However, in absence of this heterogeneity, strike-slip faulting provides an efficient mechanism for strain partitioning (compare PIV data from models A and B in Figs. 3 and 4) due to the buttress effect associated with the heterogeneity. Modeling results suggest a similarity with observations made in Iberia. The presence of lateral lithospheric variations in central Iberia may have actively contributed to the tectonic uplift of the Spanish and Portuguese Central System during the final stages of deformation (Oligocene–early Miocene), suggesting that preexisting basement faults may have played a minor role in the evolution of the E-W–trending chain.

In addition, analogue experiments demonstrate the contribution of these fault systems to the overall high topography along the Basque-Cantabrian and Iberian Ranges (Demanda Range). However, the magnitude of displacement rates varies considerably depending upon proximity to the collision border. This is in agreement with the large amount of displacement observed in the natural prototype along the Cantabrian-Ubierna corridor (>40 km) and at Somolinos (>35 km; Bergamin et al., 1996; Tavani et al., 2011).

The modeling of the vector field carried out over the surface of the analogue experiments showed dextral strike slip along the main preexisting faults, and restraining bends laterally developed at the fault tips, leading to E-W thrust faults. Moreover, several episodes of successive fault reactivation occurred during N-S shortening, involving different tectonic regimes (i.e., transcurrent to transpression). Similar conditions have been reported along fault segments in the Iberian Range, where tectonic inversion of preexisting rift-related structures has occurred (Casas-Sainz and Maestro-González, 1996; de Vicente et al., 2009). This complex scenario may explain the observed differences in the present-day topography across Iberia, providing a coherent model for strain distribution during the N-S Alpine shortening event.

Analogue experiments provide new insights into the Alpine evolution of preexisting Iberian structures developed during the formation of the Cantabrian orocline. The raising of the Cantabrian orocline between ca. 310 and 300 Ma led to important crustal thickness variations, injection of magmatic bodies, and widespread faulting along the western European Variscan belt. Analogue experiments showed such rheological variations in lithosphere strength play a key role in the nucleation of new structures. During the final stages of N-S post-Variscan convergence, the orocline was retightened, giving rise to a conjugate system of faults trending NW-SE and NE-SW to NNE-SSW, which accommodated most of the late orocline deformation. These structures strongly influenced the subsequent geological evolution of the western Iberian Massif and the eastern Mesozoic basins during the Mesozoic, as shown by the analogue experiments. The Cenozoic continental convergence between Africa and Iberia led to important geodynamic plate configuration changes that were responsible for simultaneous reactivation and tectonic inversion of previous structures, with the following important tectonic implications: (1) Major strike-slip movements ceased by the final stages of deformation (early Miocene); (2) NNE-SSW–trending structures underwent little displacement, probably absorbed by adjacent reactivation of NE-SW structures; and (3) the latter together with the NW-SE faults had a complex activity, leading to strain partitioning mechanisms that efficiently transferred the Alpine strain from the northern border of Iberia toward the plate interior. The small amount of displacement observed along some preexisting structures can be explained by the observation of movements in opposing directions or a combination of different tectonic movements. The presence of igneous bodies in central Spain may have influenced the formation of E-W–trending structures as a result of thermomechanical variations in lithospheric strength derived from different extension modes in the Iberian Massif and the eastern Iberian basins. These differences may have facilitated fault reactivation and topographic uplift. Our results contribute to clarify the role of late Variscan structures in the final evolution of intraplate mountain building during N-S Alpine convergence, providing new insights into the present-day configuration of topography in Iberia.

This work was funded by the Spanish Ministry of Science, Innovation and Universities under the project IBERCRUST (PGC2018–096534-B-100) and by Происхождение, металлогения, климатические эффекты и цикличность Крупных Изверженных Провинций (КИП; Origin, Metallogeny, Climatic Effects, and Cyclical Large Igneous Provinces; 14.Y26.31.0012; Russian Federation) to Gutiérrez-Alonso. Fruitful field discussions with M.I. Benito and P. Suárez-González are greatly appreciated. We are indebted to Daniel Pastor-Galán, the editor, and three anonymous reviewers for comments and suggestions that aimed at improving this manuscript.

Gold Open Access: This paper is published under the terms of the CC-BY-NC license.