Summary
We use a thin-shell approximation for the lithosphere to model the neotectonics of the Gulf of Cadiz, SW Iberia margin and the westernmost Mediterranean, in the eastern segment of the Azores-Gibraltar plate boundary. In relation to previous neotectonic models in the region, we utilize a better constrained structural map offshore, and the recent GPS measurements over NW Africa and Iberia have been taken into account, together with the seismic strain rate and stress data, to evaluate alternative geodynamic settings proposed for the region. We show that by assuming a relatively simple, two-plate tectonic framework, where Nubia and Eurasia converge NW-SE to WNW-ESE at a rate of 4.5-6 mm yr-1, the models correctly predict the amount of shortening and wrenching between northern Algeria-Morocco and southern Spain and between NW Morocco and SW Iberia, as estimated from both GPS data and geological constraints. The consistency between modelled and observed velocities in the vicinity of Gibraltar and NW Morocco indicates that forcing by slab sinking beneath Gibraltar is not required to reproduce current horizontal deformation in these areas. In the Gulf of Cadiz and SW Iberia, the modelling results support a diffuse Nubia-Eurasia Plate boundary, where the convergence is accommodated along NNE-SSW to NE-SW and ENE-WSW thrust faults and WNW-ESE right-lateral strike-slip faults, over an area >200 km wide, in good general agreement with the distribution of the seismic strain rate and associated faulting mechanisms. The modelling results are robust to regional uncertainties in the structure of the lithosphere and have important implications for the earthquake and tsunami hazard of Portugal, SW Spain and Morocco. We predict maximum, long-term average fault slip rates between 1-2 mm yr-1, that is, less than 50 per cent the average plate relative movement, suggesting very long return periods for high-magnitude (Mw > 8) earthquakes on individual structures.
1 Introduction and geological/geodynamic setting
The present-day boundary between Nubia and Eurasia tectonic plates changes from dextral transtensional in the Terceira Ridge (east of the Azores Triple Junction), through a dextral transcurrent (transform-type) boundary along the Gloria Fault to compressional with minor right-lateral strike-slip in northern Algeria (Fig. 1). According to recent kinematic plate models, based on GPS data, the convergence between Nubia and SW Eurasia (Iberia) is oblique, striking NW-SE to WNW-ESE (red arrows in inset of Fig. 1), and occurs at a rate of 4.5-6 mm yr-1 (Sella et al. 2002; Calais et al. 2003; McClusky et al. 2003; Fernandes et al. 2003). However, the way this convergence is being accommodated in the SW Iberia margin (comprising the southern Tagus and Horseshoe abyssal plains, the Gorringe Bank and the Tore-Madeira Rise), the Gulf of Cadiz and the westernmost Mediterranean (comprising the Alboran Sea and surrounding Betic and Rif cordilleras, often designated as the Alboran domain) is still a matter of controversy (e.g. Gutscher et al. 2002; Fadil et al. 2006; Stich et al. 2006; 2010; Zitellini et al. 2009; Vernant et al. 2010).
Topography-bathymetry map of Iberia, NW Africa and the surrounding central eastern Atlantic and western Mediterranean Sea showing the instrumental seismicity (red circles). The epicentres were extracted from the catalogues of the International Seismic Centre (ISC, between 1964-2006; ), the Instituto de Metereologia de Portugal (IM, 1970-2000; ), the Instituto Geográfico Nacional (ICN, 2002-2007; ) and the Centre Sismologique Euro-Méditerranéen (CSEM, 2007-Fev 2009; ). In the periods of superposition the IM catalogue prevailed, then the ISC. Focal mechanisms are from University of Harvard catalogue for the whole area (), complemented with a compilation from smaller magnitude earthquakes from several published sources. Between 20°W-5°W, the data sets were completed from various sources. Arrows at right bottom corner show the relative movement of Africa/Nubia with respect to Eurasia at the centre of the Gulf of Cadiz, according to different authors. Black arrows deduced from geological indicators (Argus et al. 1989; DeMets et al. 1994) and red arrows from GPS data (Sella et al. 2002; Fernandes et al. 2003; Calais et al. 2003). The white square delimits the area of Fig. 2, where the recent SW Iberia-Gulf of Cadiz structural map is discussed. Inset on the left-bottom corner shows the distinct segments where the Africa-Eurasia plate boundary has been identified (thick black lines) and where is still poorly constrained (stippled area), and the relative movement of Nubia with respect to Eurasia at distinct locations of the boundary (red arrows; after Fernandes et al. 2003). Acronyms in alphabetic order: AGFZ, Azores-Gribraltar Fracture Zone; AM, Atlas Mountains; AS, Alboran Sea; BM, Betic Mountains; GC, Gulf of Cadiz; GB, Gorringe Bank; GF, Gloria Fault; GL, Gulf of Lyon; HAP, Horseshoe Abyssal Plain; IAP, Iberia Abyssal Plain; MAR, Mid-Atlantic Ridge; Pyr., Pyrenees; RM, Rift Mountains; SAP, Seine Abyssal Plain; SBS, South Balearic Sea; SoG, Straits of Gibraltar; TAP, Tagus Abyssal Plain; TM, Tell Mountains; TMR, Tore-Madeira Rise; TR, Terceira Ridge; TTF, Tell Thrust Front; VT, Valencia Trough.
The region forms an approximately 1000-km-long, and up to 400-km-wide corridor of widespread seismicity, where the location and style of the plate boundary is still poorly constrained (stippled area in inset of Fig. 1), and where focal mechanism solutions of shallow to intermediate depth earthquakes show important variations in the regional faulting trends and rupture mechanisms (Buforn et al. 2004; Stich et al. 2006; 2010; Fig. 1) with (1) predominant strike-slip, with some reverse faulting, under NW-SE compression along the west Portuguese margin; (2) reverse and strike-slip faulting in the SW Iberia margin and Gulf of Cadiz, with maximum compression in the NW-SE direction and (3) a combination of strike-slip and normal faulting in the Alboran Sea and minor thrusting in the surrounding cordilleras, under NNE-SSW to N-S compression. Beneath the Alboran Sea, a lithospheric slab extending to depths of more than 600 km has been inferred from the seismicity and tomographic images (Buforn et al. 1991; Blanco & Spakman 1993; Seber et al. 1996; Calvert et al. 2000; Wortel & Spakman 2000).
Several geodynamic models have been proposed to explain the seismicity of the Gulf of Cadiz-Alboran region, the nature and geometry of the lithospheric slab, as well as the Oligocene-Recent evolution of Alboran domain, marked by the strong clockwise/counter-clockwise rotation of the allochthonous terranes and external zones of the Rif and Betic mountain belts, (respectively) during the Miocene (Allerton et al. 1993; Platzman et al. 1993; Lonergan & White 1997), and by the apparently coeval extension of the Alboran Sea. These include (1) backarc extension due to subduction roll-back (Royden 1993; Lonergan & White 1997; Gutscher et al. 2002); (2) extension induced by the break-off of a subducting lithospheric slab (Blanco & Spakman 1993; Carminati et al. 1998); (3) delamination of subcontinental lithosphere (Docherty & Banda 1995; Seber et al. 1996; Calvert et al. 2000; Valera et al. 2008) and (4) convective removal of a thickened lithospheric root (Platt & Vissers 1989; Platt et al. 2003).
West of the Straits of Gibraltar, the Gulf of Cadiz and west Iberia margin formed as extensional and/or pull-apart rift basins between the Late Triassic and the Early Cretaceous, in a sequence of Mesozoic extensional events (Pinheiro et al. 1996; Terrinha 1998; Alves et al. 2009). The segmentation and evolution of the basins was strongly controlled by the Variscan structural framework, in particular a set of subvertical lineaments formed during the late Variscan compression (Ribeiro et al. 1979, 1990). These lineaments strike predominantly NNE-SSW to ENE-WSW and NNW-SSE to NW-SE along the western Portuguese margin (Ribeiro et al. 1979; Murillas et al. 1990; Pinheiro et al. 1996; Alves et al. 2009), and ENE-WSW in the southern margin (Terrinha 1998; Terrinha et al. 2002; Carrilho et al. 2004; Fig. 2).
Simplified tectonic map of the SW Iberia-Gulf of Cadiz region (modified from Terrinha et al. 2009). Offshore, we used the high-resolution SWIM bathymetric compilation (Zitellini et al. 2009) complemented by GEBCO (2003). See text for description of the structures and respective references.
The present tectonics of the region is dominated by the reactivation of these late Variscan lineaments as thrust and/or strike-slip faults, which delimit prominent reliefs (e.g. the Gorringe, Marquês de Pombal and Tagus Abyssal Plain thrusts in Fig. 2) and are in places associated with fault scars and important mass wasting deposits (Gràcia et al. 2003; Zitellini et al. 2004). Recently, Duarte et al. (2009) and Zitellini et al. (2009), based on high-quality multibeam swath bathymetry (Diez et al. 2006) and multichannel seismic data, identified a new set of major tectonic lineaments striking WNW-ESE between the western Horseshoe Abyssal Plain and the eastern Gulf of Cadiz (the SWIM lineaments in Fig. 2). These cut through distinct morphological domains and show, in places, evidence of recent dextral strike-slip movement (Duarte et al. 2009; Rosas et al. 2009). According to Zitellini et al. (2009), the SWIM lineaments could represent the surface expression of a ‘precursor’ plate boundary between Nubia and Iberia and, therefore, the transition from a diffuse (Sartori et al. 1994; Hayward et al. 1999) to a discrete plate boundary setting.
In this study, we apply thin-shell lithosphere finite element techniques (Bird 1999) to model the neotectonics of the SW Iberia margin, Gulf of Cadiz and the Alboran domain. Our main goals are to improve the current understanding of how the strain is accommodated along the eastern segment of the Azores-Gibraltar plate boundary (see inset of Fig. 1) and its prolongation into the western Mediterranean, and to put some constraints on the present-day location and kinematics of the plate boundary.
Thin-shell (or thin-sheet) finite element modelling has been demonstrated to be a powerful means to investigate the neotectonics of complex plate boundary settings, characterized by marked transitions in the seismotectonic regime and lateral variations in nature and structure of the lithosphere (e.g. Bird & Kong 1994; Jiménez-Munt et al. 2001; Liu & Bird 2002; Negredo et al. 2002). The models incorporate realistic boundary conditions based on plate motions, fault networks, thermally activated non-linear rheologies for the crust and mantle, and the stresses associated with horizontal pressure gradients (due to topography and its compensation), and make predictions of strain partitioning, stress orientations, velocity fields and fault slip rates (Kong & Bird 1995; Bird 1999), which may be quantitatively compared against geophysical, geodetic and geological observations.
Previously published neotectonic models for the eastern segment of the Azores-Gibraltar plate boundary have successfully reproduce the major changes in the stress regime between the Tell mountains and the Gloria Fault, and obtained a good correlation between the highest predicted strain rates and the regions of strong seismic activity (e.g. Jiménez-Munt et al. 2001; Negredo et al. 2002; Jiménez-Munt & Negredo 2003). However, these earlier studies did not have access to the GPS velocities compiled for Iberia and NW Africa over the last decade (e.g. Stich et al. 2006; Fadil et al. 2006; Vernant et al. 2010), and assumed for the relative motion between the Nubia and Eurasia the rotation poles and angular velocities calculated Argus et al. (1989) and DeMets et al. (1990), which are based on geological indicators (i.e. magnetic lineations and transform faults). As depicted in the inset of Fig. 1, these are strongly oblique in relation to the predictions from recently derived kinematic plate models based on space-geodetic solutions.
Moreover, in contrast to previous works, which either consider the entire Azores-Gibraltar Plate boundary, extending to the mid-Atlantic ridge (Jiménez-Munt et al. 2001; Jiménez-Munt & Negredo 2003), or focus in the Iberian-Maghrebian region (Negredo et al. 2002), we attribute a greater relevance to the SW Iberia margin and Gulf of Cadiz areas. We had access to a better constrained tectonic map offshore, based on recently acquired multibeam bathymetry, backscatter data and numerous high-quality multichannel seismic profiles (Terrinha et al. 2009; Zitellini et al. 2009 and references therein). As depicted in Fig. 2, the tectonic map is dominated by NNW-SSW to NE-SW and ENE-WSW trending thrusts, and by long, WNW-ESE dextral strike-slip faults, both of which have been proposed as possible sources for destructive earthquakes and tsunamis in the region (Zitellini et al. 2001; Baptista & Miranda 2003; Gràcia et al. 2003; Terrinha et al. 2003; Ribeiro et al. 2006; Stich et al. 2007; Baptista & Miranda 2009; Zitellini et al. 2009; Cunha et al. 2010). By incorporating these structures in the neotectonic models we will be able to provide some estimates of fault slip rates, and thus evaluate their relative long-term seismic hazard. Another significant novelty of this study is that we evaluate the potential effects on the neotectonics of the area of plate boundary geometries which have been recently proposed for the SW Iberia-Gulf of Cadiz-Alboran domain region, based on the tectonics (Zitellini et al. 2009), geodetic studies and seismic data (e.g. Gutscher 2004; Stich et al. 2006).
2 Methodology
We use the thin-shell finite element program SHELLS (Kong & Bird 1995; Bird 1999) to model the neotectonics of the eastern sector of the Azores-Gibraltar Plate boundary. The thin-shell lithosphere approximation implies that only the horizontal velocity components of the momentum equation are integrated and solved across a finite element grid (FEG) of spherical triangles. Body forces arising from horizontal gradients of pressure are included in this modelling. Because the angular velocity is assumed invariant with depth, each triangular element deforms by a rigid plate rotation and a uniform strain rate. The calculation of the velocity field, and associated strain, further assumes that the lithosphere behaves as a continuous medium, even in the neighborhood of plate boundaries (Kreemer et al. 2000). This is a reasonable approximation when the modelled area is significantly larger than the thickness of the elastic/brittle lithosphere (England & McKenzie 1982).
Although the velocity model is 2-D, the momentum equation is solved accounting for the vertically integrated strength of the lithosphere. In this sense, SHELLS can be regarded as a 2.5-D finite element method (Kong & Bird 1995). At each node of the FEG, the model imports the elevation and surface heat flow and calculates, iteratively, the crustal and mantle lithosphere thicknesses under the assumption of local isostasy. To compute the temperature distribution, the steady-state vertical heat conduction equation is solved assuming a constant temperature at the base of the lithosphere, Ta. The system is considered to be isostatically balanced with a 7-km-thick mid-ocean ridge at a depth of 2.7 km. Radiogenic heat production is assumed constant within the crust (i.e. an averaged value, following Negredo et al. 2002), and negligible in the lithospheric mantle. The default density and thermal parameters used throughout this study are summarized in Table 1, and are similar to those used in previous thin-sheet neotectonic models of the region. (Negredo et al. 2002; Jiménez-Munt et al. 2003).
Model parameterization.
is the anelastic strain rate tensor, T is the absolute temperature, z is the depth and A, B, C and n are rheological parameters. The term within the square root is the second invariant of the strain rate tensor, and parameter B in the crust is given by Q/nR, where Q is the activation energy and R is the gas constant (Bird 1989). The rheological parameters are adopted from studies in the region (Negredo et al. 2002; Jiménez-Munt & Negredo 2003).3 Model construction
In SHELLS, a model is essentially defined by the following inputs: (1) the geometry of the FEG and the traces and dips of the active (or potentially active) faults; (2) the boundary conditions and (3) the topography and heat flow data, which will determine how the thicknesses of the crust and lithospheric mantle vary laterally.
For this study, the modelled region was defined between parallels 30 °N and 45°N and meridians 15°W and 3°E. The lateral boundaries were chosen so that the modelled domain encloses the area where the geometry of the Nubia-Eurasia Plate boundary is less well-constrained (Fig. 1). To the west, the boundary is a clearly marked, E-W oriented fracture zone (the Gloria Fault). To the east, the plate boundary is also linear and well-defined, following the thrusts of the Tell Mountains in northern Algeria. The southern boundary is placed to the south of the Atlas, to allow for some convergence accommodation in this area and the northern boundary is assumed to be representative of the motion of the Eurasian Plate. The model boundaries extend well beyond the SW Iberia margin, Gulf of Cadiz and west Alboran domains, where we will focus our analysis (Figs 1, 2 and 3). In this way, we avoid the modelling results in our target area being strongly affected by the imposed boundary conditions.
(a) Finite Element Grid (FEG) built for the reference Model-1, and geometry and kinematic settings of the model boundaries. The FEG is composed of continuous (thin grey) and fault elements (thick black), weaved from a regular node space grid of ~30 km. Symbols on the fault traces represent fault type: triangle-thrust fault; short line-normal fault. Faults without symbols attached are strike-slips. The coast line is shown in orange. Thick black open triangles denote model boundaries (see discussion in the text). The input, regional velocity field is depicted in red along the Nubia Plate model boundaries. Fixed model boundaries have zero velocity. Structural elements (see text) acronyms in alphabetical order: AC, Alpujarran Corridor; Al-Pa-Ca, Alhama-Palomares-Carboneras system; AR, Alboran ridge; BTF, Betics Thrust Front; CA, Cadiz-Alicante lineament; GB, Granada Basin; GF, Gloria Fault; HATF, High Atlas Thrust Front; MATF, Middle Atlas Thrust Front; RTF, Rif Thrust Front; TTF, Tell Thrust Front; YF, Yussuf fault. b) Zoom of the FEG in the SW Iberia Margin and Gulf of Cadiz. Acronyms in alphabetical order: CF, Cadiz fault; GF, Guadalquivir fault; Go-T, Gorringe thrust; HSF, Horseshoe fault; MPF, Marquês de Pombal Fault; PB-T, Portimão Bank Thrust; SM-1,2,3, SWIM lineaments 1, 2 and 3; TAP-T, Tagus Abyssal Plain Thrust.
3.1 Grid geometry and structure
The model FEG should be defined as giving a good compromise between computing efficiency and accuracy in representing the structural features in the target area. We found that an initial grid spacing of ~30 km is sufficient to reproduce the lateral variability of the lithosphere in the modelled region. In areas of structural complexity, and considering the size of the modelled faults (we only model structures with regional expression) we reduced the elements size so that locking at the fault tips and tight corners do not reduce significantly the motion along the faults, or change the patterns of strain. Fig. 3(a) shows the FEG that corresponds to our starting (or reference) model, Model-1. It comprises 4412 spherical triangles, delimited by 7936 continuous elements and 350 fault elements. The node spacing was, in places, reduced to less than 5 km, namely off SW Iberia and in the Gulf of Cadiz (Fig. 3b; see also the structural map in Fig. 2).
Fault elements are defined by their length, strike and dip. Due to the lack of geological and geophysical constraints on the gradient of most fault planes, dips of 30° and 65° were tentatively attributed to thrust and normal faults, respectively. Except in the case of vertical faults, which the model assumes as strike-slip faults, the dip assignment does not force the type of movement in the fault; that is, the fault motion sense and slip rate are model outputs. The reader is referred to the appendix in the study by Bird & Kong (1994) for a detailed explanation on the incorporation of fault elements in the thin-sheet approach, which takes into account the strengths of both the frictional and the creeping layers.
East of the Strait of Gibraltar, Model-1 is based on the previously published neotectonic models of Negredo et al. (2002) and Jiménez-Munt & Negredo (2003). At the scale of these models, a simplification of the local geology is commonly assumed, where disrupted, subparallel fault traces are often represented by long, continuous fault zones which accommodate most of the deformation. This is the case, for example, of the Cadiz-Alicante strike-slip fault and the large thrusts fronts in the Atlas mountains, in the Rif chain and the Tell Mountains (Fig. 3a). The thrust front along the Betic chain, continuous for more than 600 km, is also an exaggerated representation of several mapped compressional structures, bordering both the most external (thin-skinned) and the internal units of the mountain belt (Negredo et al. 2002).
Logically, the slip rates predicted by the models along these faults will depend not only on their extent, but also on their orientation relative to the present-day stress field and connectivity with other fault systems. In this sense, the imposed lateral continuity of the fault zones should be understood as an end-member assumption to assess their importance in the neotectonics of the region. Also considered in the southern Iberia region are the Granada basin NNW-SSE trending normal faults and the NE-SW to NNE-SSE Alhama-Palomares-Carboneras fault system (Negredo et al. 2002 and references therein). In the Alboran Sea, the Nerja and Yussuf NW-SE trending faults and the NE-SW Alboran Ridge (see Fig. 3a for location) have been mapped from wide-angle seismic data (Comas et al. 1999) and show evidence of recent tectonic activity (e.g. Watts et al. 1993; Alvarez-Marrón 1999).
Fig. 3(b) zooms on the Gulf of Cadiz and SW Iberia margin, where our models are based on recently acquired high-resolution bathymetry, backscatter and seismic data (Fig. 2). In our reference Model-1, relatively long (up to 200 km), although disrupted, SWIM segments are assumed in areas where there are some evidence of their recent and/or past activity (Duarte et al. 2009). Also included are the NE-SW thrusts off southern Portugal, linked through possible NW-SE transfer zones (Zitellini et al. 2001; Terrinha et al. 2003; Cunha et al. 2010), the Gorringe Bank thrust, the easternmost segment of the Gloria Fault (Fig. 3a), the Cadiz Fault and the Portimão Bank and Guadalquivir thrust fronts (see also structural map of Fig. 2). As discussed in Section 1, all these structures show evidence of recent (Pliocene-Quaternary) activity and are, in places, associated with seismic activity. We explicitly did not consider the Gulf of Cadiz Accretionary Wedge in the models, a prominent feature of the tectonic map presented in Fig. 2, because it consists of a sequence of stack thrusts soling along a basal, near top Cretaceous decóllement horizon, not affecting the underlying crust/lithosphere. Variations to this model setup will be discussed in Section 5.2, together with the modelling results.
3.2 Model boundaries
Another requirement of the SHELLS modelling approach is the definition of the model boundary types and forces. The models can be tectonically loaded from below, through asthenosphere coupling and/or from the lateral boundaries. In this work, we will assume that no horizontal shear traction is exerted on the base of the lithosphere and consider only the effects due to loading from the lateral boundaries, that is, the stresses resulting from the relative plate motions and their interactions.
Several kinematic plate models have been recently proposed, which place constraints on the relative motion between the Nubian plate and Iberia (e.g. Sella et al. 2002; Calais et al. 2003; McClusky et al. 2003; Fernandes et al. 2003). Fernandes et al. (2003; model DEOS2K), for example, based on GPS data predicts a convergence rate between Nubia and Europe of 4.3 mm yr-1 with a direction of N70°W in the center of the Gulf of Cadiz (9°W, 36°N). This convergence is significantly tilted in relation to the one predicted by the kinematic global plate model NUVEL-1A of DeMets et al. (1994; 4.1 mm yr-1 striking N53°W), commonly used as a reference in geodynamic studies (see inset in Fig. 2 for the velocity vectors associated with several proposed kinematic models) we did not test the recent MORVEL global plate model (DeMets et al. 2010) because the rotation pole is almost identical to the one in model NUVEL-1A.
A number of recent works, also based on space geodetic observations, confirm the DEOS2K model (e.g.Calais et al. 2003; McClusky et al. 2003; Kreemer et al. 2003) and predict the location of the Nubia-Iberia Euler pole of rotation, in average, between 10°S-1°N and 28°-15°W, that is, more than 20° to the south of the pole inferred by of DeMets et al. (1994; 21°N, -21°W; kinematic model NUVEL-1A). An important difference between the NUVEL-1A and the recently proposed GPS based models is that NUVEL-1A assumes a single African block, instead of considering the Nubia and Somalia as two independent lithospheric plates with an opening rate of up to 7 mm yr-1 along the East African Rift (Fernandes et al. 2004). One of the aims of this work is thus to evaluate the consequences of using different kinematic models on the calculated velocity fields and strain rates in the study area.
For simplicity, the Eurasia Plate will be assumed as reference (fixed plate) throughout this study. As depicted in Fig. 3(a), three different types of boundaries may be defined in the models (1) fixed, as established along the northern limit of the model and north of 37.7°N in the western limit (thick, open triangles); (2) ‘free’ (in the sense that motion is allowed), along the eastern boundary north of 36.8°N and along the western boundary between 36°N and 37.7°N, that is, in the vicinity of the Gloria Fault and (3) moving according to a particular kinematic plate model, as defined along the south, southeast and southwest boundaries of the model (red arrows). In the case of the ‘free’ boundaries, they are only subjected to lithostatic normal traction, with no shear traction. The western boundary defined in the model presented in Fig. 3 allows for a broad area (~200 km wide) of wrenching around the Azores-Gribraltar Fracture Zone (see inset of Fig. 1).
3.3 Lithospheric structure
As discussed in Section 2, the thermal and mechanical structure of the lithosphere is computed at each node of the FEG from the elevation, the surface heat flow and the rheology of the crust and mantle lithosphere, assuming local isostasy, a steady state thermal regime and the thermal-density parameters of Table 1. We describe below the utilized data sets and compare, with considerable detail, the model predictions with the available evidences and previous lithospheric models.
The topography data set used was the GEBCO 1 arc-minute grid, which provides a continuous digital terrain model over oceans and land. Offshore, the data is mainly based on digitized bathymetric contours from the GEBCO Digital Atlas, together with shallow-water soundings where available (GEBCO, British Oceanographic Data Centre 2003). Onshore, the elevations are from the GLOBE digital elevation model (GLOBE Task Team & others 1999). The data has subsequently been re-sampled at 5 arc-minute grid spacing for usage in SHELLS, as this resolution provides sufficient information for the modelling. All gridding and re-sampling operations have been performed using the open source GMT software (Wessel & Smith 1998).
For the surface heat flow, we used an updated and re-gridded version (5 arc-minute spacing) of the data set compiled by Jiménez-Munt & Negredo (2003). The utilized data set includes the measurements available from the global data set of Pollack et al. (1993), complemented by data in the Iberia Peninsula, its margins and the Western Mediterranean compiled by Fernàndez et al. (1998). A few measurements offshore Galicia (~42°N; Louden et al. 1997), the SW Portuguese margin and Gulf of Cadiz (Grevemeyer et al. 2009) have also been included in the compiled grid. In the Atlantic, offshore Iberia and Morocco, the lack of data was compensated by assuming an age-heat flow relationship from the cooling plate model of Parsons & Sclater (1977).
Figs 4(a) and (b) show the model predicted crustal and lithospheric thicknesses, respectively. In general, crustal thicknesses increase from <10 km in the distal West Iberia margin to ~30 km along the coast, in good agreement with the results from numerous seismic refraction/wide-angle experiments, both offshore (e.g. Pinheiro et al. 1992; Whitmarsh et al. 1996; González et al. 1999, 2001; Dean et al. 2000; Afilhado et al. 2008) and onshore (Córdoba et al. 1988; Banda 1988; Díaz et al. 1993; Matias 1996). A similar structure is also recovered off NW Morocco, where slightly higher crustal thicknesses are found offshore (Contrucci et al. 2004), and along the northern Spanish coast, between 5°W and 9°W approximately (Álvarez-Marrón et al. 1997 and references therein).
Model calculated lithospheric structure, assuming local isostasy and a steady state thermal regime (see Section 2 of the text for details). The main thermal-density parameters are provided in Table 1. (a) Crustal thickness, with contours every 2.5 km. (b) Total lithopsheric thickness, with contours every 10 km.
Between the Gulf of Cadiz and the Horseshoe and the Seine abyssal plains, the model shows a westward thinning of the crust from >25 to <10 km. Near the Portuguese and Spanish southern coasts, and towards the west, this crustal structure is consistent with that constrained from multichannel and wide-angle seismic profiles (Fernàndez et al. 2004; Rovere et al. 2004). Towards the central Gulf of Cadiz, however, it is likely that the model predicted crustal thicknesses include most of the thick olisthostrome body; see compilations of seismic data by Thiebot & Gutscher (2006) and Gutscher et al. (2009). In the Alboran Sea, the model predicted crustal thicknesses are consistent with those compiled by Gutscher et al. (2009) and comparable with the models inferred by Torné (2000) and Fullea et al. (2006), based on combined forward modelling of heat-flow, gravity, geoid and elevation data.
In the west and central Iberia Peninsula, as well as under the surrounding mountain ranges, the predicted crustal thicknesses are also in good general agreement with the available constraints from both seismic data (Choukroune 1992; ILIHA DSS Group 1993; Matias 1996; Pedreira et al. 2003) and tomography models (Gurría & Mezcua 2000; Souriau et al. 2008; Díaz & Gallart 2009). For example, the wide-angle/refraction seismic data show a progressive thickening of the crust between the northern Portugal and Galicia and the central Iberia System, from 30 to 34 km, approximately, and a relatively constant 30-km-thick crust in the South Portuguese and Ossa Morena terranes (Matias 1996; ILIHA DSS Group 1993). The larger discrepancies are observed in areas of large sediment accumulations, such as the high Douro and Tejo basins, to the north and south of the Spanish Central System, respectively. In both cases, however, the differences in relation to our model do not exceed 5 km (Díaz & Gallart 2009 and references therein).
In North Africa, the seismic constraints on the crustal structure are relatively scarce, or even inexistent, as in the Sahara Craton. Along the Atlas mountains, however, our model predicts a crustal root which is in average 5 km thicker than that constrained from a seismic explosion event realized along NNW-SSE, NW-SE and E-W profiles (Wigger et al. 1992). Combined 3-D topography, geoid and gravity modelling also suggest a crustal thickness varying from 34 to 38 km between the middle and the High Atlas (Fullea et al. 2007).
Fig. 4(b) shows that a regional lithospheric thickness of 85-95 km is recovered under central and west Iberia Peninsula as well as along the Atlantic margins of Iberia and Africa. Onshore, these values are slightly higher than inferred from shear wave velocity modelling (80 km; Badal et al. 1993) and deep seismic sounding investigations (85-90 km; Díaz et al. 1993). The relatively short wavelength variations observed in the map of Fig. 4(b) are either due to local crustal thickening (e.g. the volcanic seamounts and the Galicia Bank in the West Iberia margin and continental Iberia mountain ranges), important heat flow gradients (e.g. between the Gulf of Cadiz and the Alboran Sea, and across the Pyrenees), or to a combination of both (e.g. the Cantabrian mountains and Central Iberia System).
Along the continental margins, offshore and underneath Africa, previous lithospheric thickness estimates are mostly based on combined forward modelling of topography, heat flow, gravity and geoid data (e.g. Zeyen & Fernàndez 1994; Banda et al. 1995; Torné 2000; Zeyen et al. 2005; Fullea et al. 2007). In most areas, the differences between these estimates and the model presented in Fig. 4(b) are <15 km, including the Pyrenees-Catalan ranges (Zeyen & Fernàndez 1994), under the Betic chain and Alboran Sea (Torné 2000) and along the West Iberia margin (Banda et al. 1995). Large discrepancies are, nevertheless, observed in the eastern Straits of Gibraltar and Rharb Basin and in the Atlas mountains, where lithospheric thicknesses of 140-160 and 60-80 km, respectively, have been modelled (e.g. Zeyen et al. 2005; Fullea et al. 2007). The main difference is that these studies relax the heat flow constraints and are strongly influenced by geoid anomaly data. In the eastern Atlas, for example, the 2-D model put forward by Zeyen et al. (2005) overestimates the measured heat flow by up to 30 mW m-2.
The 3-D lithosphere model proposed by Fullea et al. (2007) mimics the ‘L-shaped’P-wave tomographic anomalies of Spakman & Wortel (2004) under the Gibraltar Arc, though displaced to the east by up to 100 km in the Gulf of Cadiz. On the other hand, the thick lithosphere extending from offshore NW Morocco to the Gulf of Cadiz and the Betic chain (>140 km thick) is not confirmed by recent tomography (e.g. Raykova & Panza 2010) and receiver function studies (Dündar et al. 2011). Dündar et al. (2011), in particular, based on observations from stations surrounding the Alboran Sea, central Iberia and Morocco, suggests a lithosphere-asthenosphere boundary around 90-100 km depth under NW Morocco and southern Portugal, shallowing to 60 km in the Alboran, consistent with the predictions from our model. Considering the existing observations and model predictions, the Fullea et al. (2007) thick lithosphere may be regarded as an end member model for the region.
Notwithstanding, we agree with Fullea et al. (2007) that the adopted assumption of thermal steady-state is particularly acceptable for old tectonothermal provinces. In contrast, the study area was affected by Mesozoic extension and Alpine orogeny, and therefore is likely affected by transient perturbations in the temperature distribution. Model results must therefore be interpreted as average physical conditions necessary to produce the required density distribution and/or measured heat flow, rather than the actual thermal structure. Fullea et al. (2007) suggest that the steady state assumption leads to overestimation of the actual lithospheric thickness in regions of lithospheric thinning and to its underestimation where the lithosphere is thickened.
4 Model/kinematic constraints
The thin-sheet modelling results are quantitatively compared and scored with three different sets of available kinematic and stress data: (1) stress directions; (2) seismic strain rate and (3) GPS inferred velocities. Together, these data provide a solid framework to understand the effects that different fault geometries and/or boundary conditions may have on the present-day seismotectonic regime of the region. Some changes were made in relation to previous versions of the scoring code, which improve the quality of the scoring. These include the definition of a maximum earthquake magnitude, when calculating the seismic strain rate from the recorded events, the rating of the GPS measurements according to the size of the associated error ellipses, and the definition of a subarea where the scoring is performed. For this study, we delimited the scoring area between 13°W and 0°, and between 32°N and 40°N (shaded area in Fig. 5). In this way, we avoid border effects and focus on the domains where our models are of greater relevance, namely the SW Iberia margin, the Gulf of Cadiz and the Alboran domain.
Direction of most compressive horizontal principal stress. NF, Normal faulting; NS, Normal/Strike-slip (transtension); SS, Strike-slip; TS, Thrust/Strike-slip (transpression); TF, Thrust faulting; Und., undefined fault mechanism. The grey shaded rectangle shows the area of scoring, that is, where the fit between the model predictions and observations will be measured (see text).
4.1 Stress orientations
The misfit in the principal stress directions is evaluated as the mean deviation between the observed and the computed maximum horizontal compression stress azimuths-Shmax (e.g. Jiménez-Munt et al. 2001; Jiménez-Munt et al. 2003). In total, 199 Shmax azimuths have been compiled within the scoring area (Fig. 5), of which 107 were extracted from the World Stress Map (Heidbach et al. 2008). The rest were determined from earthquake focal mechanism solutions, not included in the Heidbach et al. (2008) compilation. For scoring, the data is weighted according to the quality of the measurement, with weights varying between 1 and 5 (corresponding to qualities E to A, respectively; e.g. Jiménez-Munt et al. 2001). Over 88% of the stress data used is of class C, 4% of Class A and B and only 7% of class D.
The value of the scoring of Shmax in Neotectonic modelling was discussed in more detail by Bird (1998). He found that due to the large data dispersion in the Shmax data sets, which have a strong influence of local conditions, it is not expected that any neotectonic modelling of the lithosphere in a broad area will obtain an average error less than 25°. Moreover, keeping in mind the small range of expected variations in the values of this misfit, we can infer that even small variations of the average misfit can be significant to discriminate between different models.
Broadly, the stress indicators depicted in Fig. 5 suggest a stress regime changing from transpressive in the SW Iberia margin, to mainly thrusting in the Gulf of Cadiz, to a combination of normal fault and strike-slip in the west Alboran domain. A large dispersion is, however, observed in the data. Together with the lack of data over large areas and the local character of some measurements, not representative of the major structures in our models, high-mean misfits can be expected, exceeding 25° (Bird 1998; Jiménez-Munt et al. 2001; Negredo et al. 2002). For this reason, the quality of both the GPS and seismic strain rate scoring will be privileged in relation to that of the stress orientations.
4.2 Seismic strain
The logarithm of the seismic strain rate is compared with the logarithm of the modelled maximum (in absolute value) principal strain rate, through a normalized cross-correlation coefficient. We follow the procedure described by Jiménez-Munt et al. (2001), where at a given node of the FE grid the seismic strain rate is proportional to the sum of the scalar seismic moments of all earthquakes (Kostrov 1974), averaged with a moving Gaussian filter. The cross-correlation procedure does not compare absolute values of strain rates, but instead, it evaluates how these vary across the modelled area. This approach assumes that a constant fraction of the anelastic strain is expressed as earthquakes.
The seismic events have been compiled from four different catalogues and span between 1964 and 2007 (see caption of Fig. 1 for details on the utilization of the catalogues). Due to the moderate seismic activity in the study area, the seismic cycle for the largest earthquakes is much longer than the instrumental observation period. Because the thin-sheet methodology computes the long-term average (anelastic) strain of the models, we have limited the seismic magnitude of the earthquakes in the used catalogue to a maximum value of Mw= 6. As depicted in Fig. 6(b), this avoids the extreme localization of the seismic strain rate around a couple of locations (Fig. 6a) where very large events took place, namely the Horseshoe earthquake in 1969 (Mw= 7.8; Fukao 1973) and the El Asnam earthquake in northern Algeria in 1980 (Ms= 7.2; Ouyed et al. 1981). The cut-off value of Mw= 6, for which we assume the earthquake catalogue is complete over most of the areas investigated, is perhaps conservative, but adequate to seek a good correlation factor between the maximum modelled strain rate and the diffuse seismic strain rate.
Common logarithm of seismic strain rate calculated for a maximum earthquake magnitude (Mw) of 8.5 (a) and 6 (b). The epicentres location is shown in Fig. 1.
4.3 GPS velocities
Between 2006 and 2010, at least five studies have been published with observations from continuously recording GPS (CGPS) and survey-mode GPS (SGPS) stations in the NW Africa and SW Iberia regions (e.g. Fadil et al. 2006; Stich et al. 2006; Fernandes et al. 2007; Tahayt et al. 2008; Vernant et al. 2010). These studies greatly improve our understanding of the kinematics of this sector of the Nubia-Eurasian Plate boundary and show that (Fig. 7), first, most Iberia and Morocco form stable blocks, fixed in relation to the Eurasian and Nubian plates, respectively; and secondly, the Betic, Gibraltar and Rif systems move asymmetrically in relation to each other and the adjacent major plates.
Comparison between model calculated and GPS-inferred velocities for different boundary settings. (a) Assumes the DEOS2K plate kinematic model for the relative motion between Nubia and Iberia and a ‘free’ boundary to the west of the Iberian block. (b) As in (a) but assuming a ‘fixed’ boundary to the west in Iberia. The ‘free’ boundary segment between 36° and 37° simulates an area of distributed deformation (see text). (c) Assumes the NUVEL-1A plate kinematic model. The GPS velocities of Stich et al. (2006; orange) and Fadil et al. (2006; magenta) are residuals with respect to a stable European reference frame. The ellipses represent the error associated with the observation within a 95 per cent confidence limit. All GPS solutions are realized in the ITRF2000 global reference frame (Altamini et al. 2002).
In terms of the data distribution, the studies mentioned above can be subdivided into two groups. Stich et al. (2006) and Fernandes et al. (2007) processed mostly data from permanent GPS stations located within Iberia and southern France, whereas Fadil et al. (2006), Tahayt et al. (2008) and Vernant et al. (2010) used primarily episodic observations from Morocco. Within each group, the constrained velocity vectors may vary in their magnitude, but the directions are similar for most GPS stations. In this work, we compare the velocity fields predicted from our dynamic models to the GPS velocities constrained by Stich et al. (2006) and Fadil et al. (2006; Fig. 7), which appear as the most representative of the velocities calculated for the SW Iberia and western Morocco.
We cannot exclude the presence of coseismic and post-seismic signals in the published GPS solutions-this is observed in the solutions presented by Fadil et al. (2006) for some stations located close to the epicentre of 2004 March 6 Mw= 6.3 Al Hoceima earthquake (Stich et al. 2005). However, for most of the stations, due to the short period of observations, such contribution, if existing, is minimal and not possible to properly quantify. Therefore, we assume that the published GPS motions are not affected by co- and post-seismic signals. Because the modelled velocities are predicted long-term averages assuming a steady-state fault regime, a correction is applied to the later before scoring the models.
In essence, this correction simulates temporary (interseismic) locking of the brittle portion of the faults (Liu & Bird 2002). The misfit between modelled and the GPS constrained velocities is then expressed in terms of weighted mean error, where the weighting factor is a function of the error ellipses associated with the estimated velocities in each GPS station (Fig. 7). Because we score the models independently for each GPS data set the results will not be affected by the fact that different authors may apply distinct methodologies to compute the uncertainties in the GPS velocities.
5 Modelling results and discussion
The numerical experiments performed in this work are presented in five stages. First, we test different model boundary conditions, by comparing the computed velocity fields with the available GPS data. We then discuss the predicted sense and slip rates along the main fault systems for our reference Model-1 (Fig. 3a) in the light of a wide range of observations. Third, we test alternative tectonic settings for the SW Iberia, Gulf of Cadiz and west Alboran regions, where there are still large uncertainties on the lateral extent and connectivity between the main fault systems, their present day activity and their relation with the observed seismicity. Fourth, we score the models against the available kinematic and dynamic constraints (Shmax, seismic strain rate and GPS velocities). In Section 6, we discuss some implications for the tectonic and geodynamic setting of the region.
5.1 Boundary conditions
The model boundary conditions are essentially defined by the geometry and properties of the plate boundaries and by the kinematic model assumed to describe the relative motions between the tectonic plates. As referred in Section 3.2, we do not consider in this study the effect of basal drag, that is, which may result from mantle flow and traction at the lithosphere-asthenosphere boundary, in line with previous neotectonic modelling studies for the region (Jiménez-Munt et al. 2001; Negredo et al. 2002).
Fig. 7 shows the model predicted velocity field for three arguable boundary settings. Also plotted, for comparison, are the GPS velocities determined by Fadil et al. (2006; in magenta) and Stich et al. (2006; in orange). The models presented in Figs 7(a) and (b) both use the DEOS2K plate kinematic model of Fernandes et al. (2003) for the relative motion between Africa and Iberia (see inset of Fig. 2). They differ, however, in the type of boundary assumed west of Iberia. The model in Fig. 7(a) considers a ‘free’ boundary to the west, north of the Gloria fault, therefore allowing continuous motion in this region, similarly to that proposed by Negredo et al. (2002). These authors, however, were primarily concerned with the neotectonics of the Alboran Sea and surrounding Ibero-Maghrebian domains, and did not have access to the GPS data, only recently made available. As depicted in Fig. 7(a), this type of model predicts small relative motion between Iberia and the Maghrebian region, suggesting that Iberia is moving together with the Nubian plate, and largely overestimates the GPS velocities inferred in central and west Iberia. The number of permanent GPS sites in this region is large enough to conclude that the current velocity field is representative of the present-day kinematics of this area.
In the model of Fig. 7(b), both the northern and western Iberia boundaries are fixed to the velocity reference frame (i.e. fixed Eurasia Plate). As discussed in Section 3.2, the gap between 36°N and 38°N allows for a broad area of wrenching between Africa and Iberia, though the velocities predicted by this type of model do not change significantly if a sharp Iberia-Africa boundary is defined at the Gloria Fault (~37.2°N; Fig. 1). Despite its simplicity, this model is in general good agreement with the GPS observations and explains some relevant features such as the turn in the direction of the velocity field from NW to predominantly WNW-ESE across the Rif and Atlas mountains, and the western motion of SW Spain in relation to Europe (2.5-4.4 mm yr-1 at San Fernando according to Stich et al. 2006 and Fernandes et al. 2007, respectively), and of NW Morocco in relation to stable Nubia (~1-2 mm yr-1; Fadil et al. 2006; Tahayt et al. 2008).
Finally, we tested the NUVEL-1A plate kinematic model, which can be considered as an end-member boundary condition, that is, one which presents the greatest obliquity in relation to the recently GPS-derived plate kinematic models (see inset of Fig. 1). As depicted in Fig. 7(c), the velocities calculated with this model show systematic deviations in relation to the GPS-derived velocities in some areas>. For example, in central Morocco the model predicts the anticlockwise rotation of the velocity field to occur further to the west, and along the NW Morocco and SW Iberia coasts it underestimates the western component of the velocity field measured in several GPS stations. It appears, therefore, that a fixed western Iberia boundary, combined with a plate kinematic model which predicts a WNW-ESE movement of Africa in relation to Iberia, strongly tilted in relation to the NUVEL-1A global solution of DeMets et al. (1994), are the most suitable boundary settings for the modelling (Fig. 7b).
5.2 Reference model predictions versus observations
Fig. 8 shows the long-term average fault slip rates predicted for reference Model-1, assuming the boundary conditions presented in Fig. 7(b) and the default thermal-mechanical parameters of Table 1. Ideally, the predicted slip rates should be quantitatively compared with measurements along the fault planes. Such measurements, however, are too scarce within the studied area to obtain a meaningful scoring. Consequently, the modelling results will be mostly discussed in terms of the GPS inferred deformation, the distribution of the seismic strain rate, geological constraints when available, and comparisons with previous neotectonic models in the region.
Broadly, the slip directions in the east Alboran Sea, SE Iberia and along the Tell mountains, are consistent with those predicted in previous works and with available geological, seismologic and GPS indicators (Negredo et al. 2002; Jiménez-Munt et al. 2003; Stich et al. 2006 and references therein). The model predicts the highest shortening rates in the Tell Mountains (between 2.6 and 4.3 mm yr-1). Possibly, these estimates should be considered as an upper limit, due to the exaggerated lateral continuity of the thrusts assumed in the model (Section 3.1). Notwithstanding, they are within the range of values predicted by Dewey et al. (1989) for the last 9 myr (5 mm yr-1), and inferred by Meghraoui et al. (1996) for the Quaternary (2.3 mm yr-1). More recently, Stich et al. (2006) based on GPS observations estimated more than 3 mm yr-1 of African-Iberia convergence accommodated in the north Algerian margin. This suggests that, in the long term, most of the shortening in the region is being accommodated by faulting, consistent with the occurrence of frequent, reverse faulting, moderate to large earthquakes (e.g. Buforn et al. 2004; Stich et al. 2006). Between SE Iberia and Nubia, the model predicts up to 2.8 mm yr-1 of dextral motion, accommodated along the Yusuf and Alpujarra faults corridors (added 1.2 and 1.6 mm yr-1, respectively). This amount of wrenching is also comparable with that inferred from the GPS data by Stich et al. 2006 (~3.5 mm yr-1). In the Cadiz-Alicante fault zone, which has a less favourable WSW-ENE direction, the model predicts a much lower slip rate (~0.1 mm yr-1), except where it intersects the Alpujarra Corridor (see also Fig. 3 for faults nomenclature).
Further to the west, significant shortening is also predicted along the faults that border the Alboran Ridge (up to 2.8 mm yr-1). In this area, however, the focal mechanism solutions of shallow earthquakes indicate predominantly sinistral strike-slip and minor normal fault motion (Buforn et al. 2004; Stich et al. 2006). In fact, only one intermediate depth (~80 km depth), reverse fault earthquake has been reported in the area (Mw= 3.9; Buforn et al. 2004). The model presented in Fig. 8 is also inconsistent with the high seismicity that characterizes the Murcia and Alhoceima regions (in the Betics and Rif ranges, respectively; Fig. 1).
In the western High Atlas, the model predicts less shortening (0.1-0.8 mm yr-1) than estimated from GPS data, which varies between maximum =1 mm yr-1 (Fadil et al. 2006; Vernant et al. 2010) and 2 mm yr-1 (Stich et al. 2006). Similar amounts of shortening (1-2 mm yr-1) have also been inferred from restored cross sections in the Middle and High Atlas (Beauchamp et al. 1999; Gómez et al. 2000; Teixell et al. 2003). A plausible interpretation is that part of the shortening is being accommodating by folding and/or by crustal-scale (or lithospheric-scale) buckling. In fact, the area is characterized by moderate to low seismicity, and only one strong event (Mb= 5.9), reverse fault earthquake has been recorded (Alami et al. 2004). In the Middle Atlas, the predicted low-slip rates (<0.2 mm yr-1) are consistent with those estimated in quaternary fault scarps (0.05-0.3 mm yr-1; Rigby 2008).
Along the western Betics and Rif mountains, our reference Model-1 shows significant differences in relation to previously published neotectonic studies, where the chains have been explicitly modelled as a single orogenic belt, connected through a major thrust front arcuated under the eastern Gulf of Cadiz (i.e. an active Gibraltar Arc; Negredo et al. 2002; Jiménez-Munt et al. 2003). Such studies predict relatively high-convergence rates in the Betics and the Rif fronts (1-2 mm yr-1 approximately), combined with moderate normal faulting and dextral strike-slip (0.3-0.5 mm yr-1) west of the Straits of Gibraltar. Contrastingly, the model presented in Fig. 8 suggests low-convergence rates in the Betics and a combination of strike-slip and thrusting in the Rif mountains, in better agreement with the shallow crustal deformation inferred from focal mechanism solutions in NW Morocco (Meghraoui et al. 1996; Buforn et al. 2004).
In the Gulf of Cadiz, Horseshoe Abyssal Plain and SW Portuguese margin, our reference fault Model-1 comprises the three main types of mapped structures: (1) large WNW-ESE trending strike-slip faults (the eastern Gloria Fault segment and the disrupted SWIM lineaments, ca. 200 km long segments); (2) the Portimão Bank and Guadalquivir ENE-WSW trending, dipping to the north thrust faults and (3) the NNE-SSW to NE-SW, dipping to the east thrusts (e.g. the Horseshoe, the Marquês de Pombal and Tagus Abyssal Plain thrusts), connected along NW-SE trending, sinistral transfer faults. As depicted in Fig. 8, the model predicts a distribution of the deformation along the mapped structures, with: up to 1 and 1.4 mm yr-1 dextral slip in the northern SWIM and Gloria fault lineaments, respectively; up 0.5 mm yr-1 in the Guadalquivir thrust and a maximum 0.7 mm yr-1 in NE-SW trending thrusts off SW Portugal, combined with 0.3-0.5 mm yr-1 in the sinistral transfer faults.
The predicted strain partitioning between the different tectonic structures in the region agrees with the large dispersion of epicentres, as well as the deformation of the Pliocene-Recent sediments observed in numerous compressional and transpressional structures between south of the Horseshoe Abyssal Plain and the SW Portuguese margin (Masson et al. 1994; Hayward et al. 1999; Terrinha et al. 2003; Zitellini et al. 2004; Terrinha et al. 2009). Moreover, the sense of movement predicted in the modelled faults is consistent with the predominant stress regime inferred from intermediate to large earthquakes, which comprise shallow to intermediate depth events associated with thrust and/or dextral strike-slip faults (Fukao 1973; Buforn et al. 1988; Stich et al. 2007).
However, according to Stich et al. (2006), who projected the GPS velocities along the direction of the principal seismic stresses, the amount of shortening between NW Morocco and southern Portugal is ~1.4 mm yr-1, that is, about three times larger than predicted in Model-1 (<0.5 mm yr-1 in the Portimão Bank and Guadalquivir thrusts; Fig. 8). This suggests that, either part of the strain is being accommodated elastically within the lithosphere or, as appears to be supported by the seismicity of the region, the model underestimates the amount of thrusting in the northern Gulf of Cadiz. We also note that the distribution of the strain is not entirely consistent with that computed from the instrumental seismicity, which in the Gulf of Cadiz focus along a WSW-ENE, ~80-km-wide band to the north (Fig. 6b, for Mw= 6; see also Fig. S1d in the online Supporting Information for comparison>).
Along the SW Portuguese margin, Model-1 is consistent with the shortening estimated by Stich et al. (2006) of ~0.5 mm yr-1. Some geological constraints on the fault slip rates are also accessible off SW Portugal, where extensive and detailed seismostratigraphic work has been done in the last decade (e.g. Roque 2008). In the Marquês de Pombal Fault, for example, it appears that up to 1 km of vertical displacement accumulated within the last 5-10 myr, resulting in the condensed Tortonian (possibly base of Pliocene) to Recent sediment sequences in the fault footwall (Roque 2008). For an average fault plane dip of 24° (Zitellini et al. 2001) this corresponds to a net slip rate of 0.25-0.5 mm yr-1, that is, slightly less than the predictions from our model (up to 0.7 mm yr-1; Fig. 8). However, since the convergence between Nubia and Iberia rotated from NNW-SSE in the Middle Miocene, to WNW-ESE in recent times (e.g. Ribeiro et al. 1996), favourable to the orientation of the NNW oriented Marquês de Pombal thrust fault, it is likely that the observed displacement accelerated in recent times, even exceeding the model predictions.
In summary, our reference tectonic Model-1 is broadly consistent with the instrumental seismicity and the kinematics of the Iberia and Nubia plates inferred from GPS data in the east Alboran, Tell Mountains and west of southern Portugal. Some improvements in relation to previously proposed neotectonic models are also observed along the Rif Mountains (NE Morocco), where the model predicts a combination of reverse faulting and dextral strike-slip and in the Atlas mountains. In the central-west Alboran, Gulf of Cadiz and Horseshoe abyssal plain, however, the model is not fully consistent with the distribution of the seismicity and, in places, with the amount of shortening inferred from GPS data. In the next two sections, we will discuss the results from two alternative fault models and quantitatively evaluate the models predictions against the available kinematic and dynamic constraints.
5.3 Tectonic setting of the west Alboran, Gulf of Cadiz and SW Iberia margin
In this section, we compare the results from three distinct structural models, representative of a large series of numerical experiments that translate some of the geodynamic models proposed for the Gulf of Cadiz-SW Iberia region. Fault Model-2 (Fig. 9a) is based on the recently published work of Zitellini et al. (2009), which identified a narrow (<100 km wide) and long (>600 km) band of strike-slip deformation (the SWIM lineaments; Fig. 2), probably linking the Gloria Fault and the Rif-Tell cordilleras. As discussed in Section 1, the authors interpreted these lineaments as a precursor of a new transcurrent plate boundary between Iberia and Africa. Model-3 (Fig. 9b), on the other hand, considers that the deformation along the SWIM lineaments is restricted to the upper crust and sedimentary cover (i.e. has limited expression at the lithospheric scale), and that the deformation is primarily accommodated along seismically active structures. To build this model, we consider a major structure comprising the Portimão Bank and Gualdalquivir thrusts off south Portugal, and a NNE-SSW trending major shear zone across the Alboran; the Trans-Alboran Shear Zone (TASZ) as proposed by Stich et al. (2006). Model-3 also tests a highly segmented thrust front along the middle-high Atlas and the central-northern Betics, where the seismic strain rate is relatively low (Fig. 6).
Long-term average fault slip rates predicted for two alternative tectonic model scenarios. (a) Fault Model-2 is based on the work of Zitellini et al. (2009) who consider the development of a new plate boundary along the WNW-ESE trending SWIM lineaments. (b) Fault Model-3 is based on the mapped structures and the distribution of the seismic strain (see text). Boundary conditions applied as in Fig. 7(b). See Fig. 3 for the names of the structures.
Fig. 9(a) shows that the model which includes a long WNW-ESE lineament, continuous for more than 600 km as proposed by Zitellini et al. (2009), concentrates almost all the deformation along a narrow zone in the vicinity of the lineament. The predicted slip rates vary along this lineament from ~2 mm yr-1, east of the Horseshoe fault, to almost 4 mm yr-1 in the Horseshoe abyssal plain. North of the wrenching zone, however, the slip rates along the NE-SW thrusts that border the southwest Portuguese coast are reduced to half, in relation to Model-1, or become insignificant, as along the north Gorringe Bank thrust. Model-2 also predicts an increase in the fault slip rates across the Rif Mountains, in the continuation of the SWIM lineaments, with a stronger strike-slip component. A major difficulty of this model is, therefore, to conciliate the existence of a developing plate boundary along the major SWIM lineaments with the seismicity of the area, which is displaced 50-100 km to the north (Figs 6 S1e), and the widespread evidences of recent re-activation of thrust faults to the north of the Horseshoe Abyssal Plain.
The predicted fault slip rates in our model Model-3 (Fig. 9b) also show important differences in relation to those of our reference Model-1 (Fig. 8). The elimination of the SWIM lineaments, for example, results in an increase in the activity of the NE-SW thrusts off southwest Portugal (up to 100 per cent in the Horseshoe fault) and along the Gloria Fault. Significantly higher fault slip rates are also predicted in the Portimão Bank and Gualdalquivir thrusts (>300 per cent in relation to Model-1). Although connected, the type of faulting predicted by the model along these two structures varies as a function of their orientation, with predominantly thrusting in the Gualdalquivir and dextral strike-slip along the more E-W oriented Portimão Bank Fault. In fact, the interpretation of multichannel seismic profiles suggests that during the Betic orogeny, in the Miocene, both structures acted as dipping to the north thrusts (Gràcia et al. 2003; Terrinha et al. 2009), whereas in recent times the deformation is masked by the northward emplacement of the olisthostrome, the overlying sediments (Camerlenghi & Pini 2009) and wrenching (Terrinha et al. 2009). The mixed strike-slip and thrust regime predicted for the area is, nevertheless, in good general agreement with the seismicity of the region, characterized by strike-slip and thrust focal mechanisms solutions (Buforn et al. 2004; Stich et al. 2006).
In the Rif Mountains, Model-3 predicts similar fault slip rates to Models-1 and 2, but with a strong predominance of thrust faulting. Although inconsistent with the lack of noteworthy thrust faulting rupture events in the area (Meghraoui et al. 1996; Buforn et al. 2004), this model should not be rejected, since the time span of the instrumental seismicity may be shorter than the recurrence period of moderate-large size earthquakes. Moreover, the model predicted shortening is consistent with that inferred by Meghraoui et al. (1996) during the Pliocene-Quaternary, based on geological and geophysical (seismic and magnetic anomalies) data, of 1-2.3 mm yr-1. In the Atlas mountains, the amount of shortening predicted by this model, which assumes a highly disrupted thrust front, is significantly less than in Models-1 and 2 (Figs 8 and 9a, respectively) and in apparent disagreement with the amount of shortening predicted from geological constraints and GPS velocities.
As mentioned above, Model-3 also simulates the TASZ as put forward by Stich et al. (2006), linking the Palomares and Carboneras faults in SE Spain to the N’kor fault in northern Morocco, across the Alboran ridge. Although the continuity between the mapped structures is a model assumption, it is in general agreement with the left lateral strike-slip motion observed in some of the main structures, such as the Palomares and Carboneras faults in SE Spain (Martínez-Díaz & Hernández-Enrile 2004; Montenat et al. 1990) and the N’kor Fault in north Morocco (Meghraoui et al. 1996), and explains a number of features in the seismicity (Stich et al. 2006,2010 and references therein), including (1) the dominance of left-lateral strike-slip events across the TASZ; (2) occasional normal and thrust faulting events, both in SE Spain and along the Alboran Ridge, enhancing the importance of fault interaction and local stress transfer and (3) the occurrence of mainly reverse fault to strike-slip events in the easternmost sector of the Cadiz-Alicante fault zone. The model predicted maximum strike-slip rates, of ~2 mm yr-1along the TASZ, between the Alpujarra corridor and the Alboran ridge, are also consistent with that estimated from GPS data (Stich et al. 2006).
Finally, Model-3 considers the Lower Tagus Valley and Nazaré faults. The Lower Tagus fault system has been invoked as responsible for moderate, but continued seismic activity (Fig. 1), and some important historical earthquakes, namely the 1531 (estimated Ms= 7.1; Sousa et al. 1992) and the 1909 (estimated Ms= 5.9; Teves-Costa et al. 1999) events, which destroyed a large part of Lisbon and Benavente, respectively. Vilanova et al. (2003) also suggested that this fault might have been re-activated during the 1755 Great Lisbon Earthquake, which had its source area offshore SW Portugal. Cabral (1995), on the other hand, based on geological criteria, estimated the Pliocene-to-Recent activity of the Lower Vale do Tejo fault in 0.05 and 0.1 mm yr-1, that is, consistent with the model predictions.
5.4 Models scoring
The results discussed above show clearly that the models predicted fault slip rates and the style of faulting depend on the assumed structural framework as well as on the faults length and segmentation. This has important implications for the geodynamic setting, the stress regime and the seismic hazard of the region. A way to quantitatively test the consistency of the different models with observations is, therefore, critical. In this study, we use three different types of constraints to score the models (see Section 4 for discussion of the data sets): (1) the seismic strain rate; (2) the direction of the maximum horizontal compression (Shmax) and (3) the GPS velocity data.
5.4.1 Seismic strain rate and Shmax
Fig. 10 compares the scoring of the three fault models discussed in the previous sections in terms of the seismic strain rate and Shmax. For each model, the fault friction coefficient (fc) is varied systematically between 0.01 and 0.85. The results show little sensitivity for fc > 0.4. For lower values, the misfit between computed and observed stress direction progressively increases, whereas the correlation with the seismic strain rate is improved for 0.03<fc<0.2, and then, in Models 1 and 2, slightly decreases for 0.01<fc<0.03. As discussed in Section 4.2, Shmax observations show high dispersion and are scarce over large areas, particularly in the deep offshore (Fig. 5). For this reason, we privilege the improvement in the seismic correlation and infer that for the modelled region 0.05<fc<0.08. This interval of fc values lies between the optimal friction coefficients found by Negredo et al. (2002) for the Ibero-Maghrebian region (0.04-0.06) and those obtained by Jiménez-Munt et al. (2001) and Jiménez-Munt & Negredo (2003) for the Azores-Gibraltar system (0.1-0.15). As suggested by Negredo et al. (2002), the low values of fc may be indicative of a lower strength of faults in continental domain.
Scoring of the three tested tectonic models, for a varying fault friction coefficient (fc), against seismic strain (right axis) and stress direction data (left axis). The improvement of the fit between modelled and observed is always upwards in the graphic. Model-3 provides the best fit to both the seismic and stress orientation data. Grey bar highlights best interval of fc values (see text).
As expected, the scoring results in Fig. 10 show that Model-3 correlates better with the computed seismic strain rate than Models 1 and 2. It is worth noticing, nevertheless, that a significant improvement, of approximately 14 per cent, is achieved in relation to our reference Model-1 (based on mapped elements), and of 40 per cent in relation to previously published neotectonic models for the region, which either considered a more restricted (Negredo et al. 2002) or a much larger (Jiménez-Munt et al. 2001; Jiménez-Munt & Negredo 2003) area for the scoring. The distribution of the model computed strain depicted in Fig. 11(a) shows, in fact, a remarkable similarity with the observed (for Mw= 6; Fig. 6b). The main discrepancies are observed within two relatively restricted areas, around Lisbon and in the central Alboran, where there are no instrumental earthquakes to justify the model predictions. We recall that the seismic coupling is assumed uniform all over the model, and that this assumption may fail in the Alboran domain, where a very thin lithosphere is inferred (<50 km; Fig. 4).
(a) Logarithm of smoothed strain rate predicted by fault Model-3, assuming a fault friction coefficient of 0.05. Arrows are predicted velocities with respect to Eurasia. (b) Shmax orientations and stress regime predicted by fault Model-3. Areas A to F (excluding E) are highlighted to compare Model-3 predicted stress azimuths with those inferred by Stich et al. (2006) based on seismic moment tensors. Computed seismic strain and observed principal stress orientations, for comparison, are depicted in Figs 6 and 5, respectively. In Fig. S1 (Supporting Information), we show the smoothed strain rate and the predicted principal stress directions for Models 1, 2 and 3.
The scoring in terms of Shmax also favours Model-3. The improvement in relation to our reference Model-1 is mainly noticed in the northern Morocco-Alboran region and, particularly, in the eastern Horseshoe Abyssal Plain-SW Iberia (between parallels 35°N and 37°N parallel and meridians 10°W and 11°W), where Model-3 is consistent with a more northerly oriented observed Shmax (Figs 5 and 11b; Fig. S1h-i).
Despite this, relatively large errors in the computed mean stress azimuth are obtained in all the models, similar to those obtained in previous neotectonic studies (Negredo et al. 2002; Jiménez-Munt & Negredo 2003). Large departures between the observations (Fig. 5) and the predictions from our best fit Model-3 (Fig. 11b), both in terms of the direction and predicted stress regime, are noticed, for example, along the southern Spanish and Portuguese coasts. As discussed in Section 4, however, the great majority of the available stress directions are computed from seismic events, often of small magnitude, which may reflect local interference patterns (e.g. Stich et al. 2010). For this reason, we have also compared our Model-3 predicted stress azimuths with those computed by Stich et al. (2006), based on moment tensor mechanisms inferred from earthquakes of Mw= 3.5. Four areas have been considered in the analysis, which are shown in Fig. 11(b). Area E of Stich et al. (2006), which comprises the Tell Mountain in northern Algeria, was not considered because the stress tensor is only poorly constrained in their work.
For the remaining areas, we obtained the following average deviations between Stich et al. (2006) and Model-3 Shmax (computed assuming Fisher statistics): A-17.1°; B-16.9°; C/D-24.9°; F-2°. Except for area A, all the other deviations are in the counter-clockwise direction. This represents an improvement of ~45 per cent in areas A and B, and of 17 per cent in areas C-D, in relation to the previous scoring, where inferior quality stress measurements were dominant. Notwithstanding, all of our models fail to reproduce the predominantly extensional stress regime in area B (Fig. S1i), where a large proportion of the computed stress tensors are associated with normal faulting (Stich et al. 2006, 2010).
Fig. 10 also shows that Model-2 has the poorest seismic correlation coefficient and the greatest error in the direction of Shmax. This is an important result, which argues against the interpretation of the mapped SWIM lineaments as mature lithospheric-scale features, extending continuously between Gloria fault and the eastern Gulf of Cadiz; in which case they would accommodate most of deformation associated with the oblique Iberia-Nubia convergence along a narrow band (~50 km wide; Fig. S1e). In fact, the progressive improvement in the scoring between Model-2 and Model-1, and then between Model-1 and Model-3, suggests a relatively minor role for the SWIM structures in the present day geodynamic context of the region.
5.4.2 Modelled velocities and GPS data
One of the novelties of this study is the usage of the recently released GPS data to constrain the modelling results. In Fig. 12, we compare the calculated velocity field for tectonic Model-3 with the GPS velocities derived by Stich et al. (2006) and Fadil et al. (2006; Fig. 12a), and analyse the mean deviations between observed (GPS-derived) and predicted velocities for the three fault models considered here, both in terms of the absolute value (Fig. 12c) and the vector azimuth (Fig. 12d). The GPS scoring is consistent with the one obtained for the seismic strain rate and Shmax, in that: (1) 0.5<fc<0.8 provides the best overall adjustment to the data, if all the models are considered; (2) Model-3 is the one that best explains the GPS observations, with a significant improvement in relation to the velocities determined by Stich et al. (2006; circles in Fig. 12c and (3) Model-2 is associated with the largest computed errors.
Predicted model velocities and comparison with GPS data. (a) Comparison between the predicted Model-3 velocities (grey arrows) and GPS velocities calculated at a number of stations in NW Africa and Iberia, according to Fadil et al. (2006; magenta arrows) and Stich et al. (2006; orange arrows). Ellipses show 95 per cent confidence limits. The thick black line is a profile along which we plot the GPS and model velocities for comparison in (e) and (f). In (b) we overlay the structural framework of Model-3 on the calculated velocity, to discuss the competing effects of the faults and lithospheric structure (see text). (c) and (d) show the mean deviation and azimuth error, respectively, between the models predicted velocities and the GPS data, as a function of fc. Tectonic models 1, 2 and 3 are shown in Figs 8 and 9. Model-3 (red) is clearly associated with the lowest mean deviation from the GPS data. The dashed lines show the scoring for Model-3 when the GPS stations in central and eastern Iberia are excluded (see text). (e) and (f) show the parallel and normal horizontal components, respectively, of both the GPS (Stich et al. 2006) and Model-3 velocities (red lines), projected along a NW-SE profile (thick black line in a). At each GPS station the one standard deviation error bar is shown. Dashed black lines give the relative motion of Nubia calculated along the profile. The dashed red line in (f) is shifted of -0.7 mm (see text for further explanation).
From the map in Fig. 12(a), we observe that the velocities calculated with our best-fit Model-3 are in general consistent with the GPS data, particularly in the structurally complex region of northern Morocco and SW Iberia, where the model reproduces both the progressive deceleration towards the north and the greater E-W component of the velocity field around the Straits of Gibraltar and northernmost Morocco. Moreover, the large azimuth errors (>20°) between the models predicted velocities and those calculated by Stich et al. (2006; red circles in Fig. 12c) can be attributed to the large deviations observed in central and eastern Iberia, where the modulus of the estimated GPS velocities are, in fact, small when compared with the obtained error ellipses (i.e. they are below the computed accuracy for the stations and can be considered as noise residuals). When these GPS station are ignored, the computed mean azimuth error decreases to about a third, and an improvement on the order of 12 per cent is achieved in the mean deviation value (dashed red lines in Figs 12d and c, respectively).
We further evaluate the quality of the fit between Model-3 predictions and the velocities estimated at the permanent GPS stations in northern Morocco and SW Iberia along a NW-SE profile, close to the direction of maximum principal stress tensor (thick black line in Fig. 12a). The profile follows as close as possible the observation sites, thus minimizing the errors in the projection of the observations to the profile. In Figs 12(e) and (f), we compare the components of the modelled (thick red line) and inferred GPS velocities parallel and normal to the profile, respectively. Fig. 12(e) shows a particular good fit to the velocities parallel to the profile. These imply ~1.8 mm yr-1 of shortening between northern Morocco and stable Nubia and ~1.4 mm yr-1 of shortening between northern Morocco and southernmost Portugal, thus supporting the thrusting rates inferred by Model-3 along the Rif and Guadalquivir-Portimão Bank fault systems (Fig. 9b). It is also worth noticing that Model-3 correctly predicts the San-Fernando velocity in SW Spain, with an intermediate velocity between Eurasia and Nubia (Fig. 12a).
Normal to the profile (Fig. 12f), our model (solid red line) shows a systematic shift in relation to the observations of ~0.7 mm yr-1, which we have corrected for (dashed red line). A similar correction has also been applied to the Nubia referential for consistency (dashed black line). After the shift has been applied, large deviations from the model are still noticed in relation those GPS stations located farther away from the profile, though within the one standard deviation error. Despite this, Fig. 12(f) shows that Model-3 correctly predicts the total variation between northernmost Morocco (Tétouan, TETN-GPS station in Fig. 12) and west of Lisbon (Cascais, CASC-GPS station), which corresponds to ~2.1 mm yr-1 of right-lateral motion, if we take into account the variation of the relative normal motion of Nubia along the profile (~0.4 mm yr-1 between the locations of Tétouan and Cascais). This amount of wrenching is about twice that predicted by Stich et al. (2006; who did not consider the variation of Nubia relative motion along the profile), but consistent with the slip rates predicted in both our Models-1 and 3 (though along distinct fault structures; Figs 8 and 9(b), respectively).
Fig. 12(b) shows the calculated velocity field overlain on our best field structural model and estimated fault slip rates. Comparisons with the velocity field calculated for our reference Model-1 (Fig. 7b) and Model-2 (Figs S1j and k) illustrate the competing effects of the tectonic framework and lithospheric structure on the velocity field. For example, across the TASZ, which in Model-3 delimits the Alboran Sea to the east, we observe an anticlockwise rotation of ~30° in the velocity field, associated with prominent left-lateral shear. As discussed above, the inclusion of the TASZ in Model-3 concords with the distribution of the seismic strain rate and the focal mechanism solutions between southern Spain and northern Morocco. An anticlockwise rotation of the velocity field is also predicted in Models-1 and 2, where the TASZ was not considered, but occurs 50-100 km further to the west and over a wider area, controlled by the contrasting lithospheric structure (thus strength) across the narrow Alboran Sea (Fig. 4) and by the faults that delimit the Alboran domain to the north and south (Fig. 3).
Other structures which appear to be associated with rotations in the velocity field include the Tell Mountains thrusts, the Yussuf corridor east of the Alboran Sea (in this case a mild clockwise rotation), the ENE-WSW to NNE-SSW thrusts off SW Iberia, and the Atlas thrusts (Fig. 12b). In the later, the anticlockwise rotation is stronger in Models-1 and 2 (Figs 7 and S1j-m), in better agreement with the GPS data. Fig. S1(k) also shows that in the case of Model-2, which considers the continuous SWIM strike-slip alignments between the Gulf of Cadiz and the Gloria Fault, a strong attenuation of the velocity field is predicted between 35°N and 36°N. Unfortunately, there is no available direct measurement offshore to confirm the models predictions, but Model-2 velocities around the Straits of Gibraltar and SW Spain are significantly smaller (up to 50 per cent) than those inferred from GPS data (Fig. S1). This is clearly reflected on the poor Model-2 scoring against the GPS constraints (Figs 12c and d).
6 Tectonic and geodynamic implications
The modelling results discussed in Section 5 show that a seismotectonic framework which assumes a diffuse Nubia-Eurasia Plate boundary between the Gloria fault and northern Algeria explains the existent kinematic and dynamic constraints, and is consistent with the amount of shortening and wrenching estimated between northern Morocco-Algeria and southern Spain-Portugal, the WNW motion of the west Alboran domain and SW Iberia in relation to Eurasia, and the left-lateral shearing along the TASZ, nearly perpendicular to the regional plate convergence.
The modelling kinematic assumptions are straightforward, as we only input the relative motion of Nubia in relation to a fixed Eurasia plate, as derived from recent geodetic data (Calais et al. 2003; McClusky et al. 2003; Fernandes et al. 2003). Therefore, although the modelling results are not irreconcilable with ongoing subcrustal regional processes, they suggest these may have a minor role in explaining, for example, the differential westward motion of the Alboran domain. As depicted in Fig. 12(a), the deviations between modelled and GPS-derived velocities in the vicinity of Gibraltar (including the GPS stations of Tetouan, Ceuta and San Fernando), as well as in NW coast of Morocco, are <1 mm yr-1. This argues against the hypothesis of a tectonic regime dominated by a presently active, eastward dipping subduction beneath Gibraltar, with slab roll-back and hydrostatic pull of the overlying plate, as put forward by Gutscher et al. (2002).
We have also considered the possibility of an active Gibraltar Arc, linking the Betic and Rif thrust fronts and explicitly separating the Alboran and Atlantic domains (Model-4 in Fig. S2b, Supporting Information). Similar to that predicted by Negredo et al. (2002), this tectonic setting results in a strong localization of the strain along the Betic chain thrust front (Fig. S2d), in disagreement with the distribution of the seismic strain rate (Fig. 6). In relation to tectonic Model-3, the scoring calculated for this model (assuming fc= 0.05; Fig. 10) is lower by approximately 15 per cent. A weaker correlation is also obtained between the model predicted velocity field and the GPS observations, by 22 and 6 per cent in relation to the Stich et al. (2006) and Fadil et al. (2006) data sets, respectively. The larger departures are observed along southern Spain and around the Straits of Gibraltar, where the model predicted velocities fail to reproduce the dominant E-W direction of the GPS data (Fig. S2d). This model also predicts much lower thrust rates along the Rif Mountains, and significantly higher (almost twice the amplitude) slip rates along the Guadalquivir-Portimão faults (Fig. S2b), thus largely underestimating the amount of shortening within the Nubia plate inferred from GPS data, and overestimating it between Morocco and southern Portugal (Stich et al. 2006; Section 5.4.2).
A feature of the GPS data which the models fail to reproduce, and that may, in fact, be indicative of coupled subcrustal (or sublithopsheric) processes and crustal deformation, is the SW motion of the Rif Mountains in relation to the Nubia Plate, and of sites in the Betics in relation to Eurasia, that is, normal to the relative motion between Nubia and Eurasia (Fadil et al. 2006; Stich et al. 2006; Vernant et al. 2010). The main reason why we did not attempt to adjust the models to these observations is that they are localized and possibly controlled by remnant lithospheric processes, associated with the NW-SE Africa-Eurasia convergence during the Eocene-Miocene. Platt & Houseman (2003), for example, argue that the tomographic data and intermediate to deep seismicity in the Alboran domain are consistent with a laterally propagating Rayleigh-Taylor instability around a slab of delaminated subcontinental lithospheric, whose downwelling initiated ca. 27 Ma (Platt et al. 1998). Pérouse et al. (2010), on the other hand, used dynamical modelling techniques and showed that the observed southward motion of the Rif and Betics can be reproduced by applying a small (100 × 50 km) traction patch to the base of an elastic plate in the external Rif zone. Based on the size, location and orientation of the patch the authors associated the traction with delamination and roll-back of subducted African lithospheric mantle.
In the Gulf of Cadiz and SW Iberia, our results evidence the inconsistency of the kinematic data with a model dominated by WNW-ESE trending long lineaments (up to 600 km), extending between south of the Gorringe and the eastern Gulf of Cadiz (the SWIM lineaments, Duarte et al. 2009; Zitellini et al. 2009). The great spatial continuity of these lineaments has been suggested based on morphological criteria (high-resolution multibeam bathymetry) and Zitellini et al. (2009) argued they represent an emergent plate boundary between Nubia and Iberia, linking the Gloria Fault to the Rif-Tell Mountains.
Such a tectonic model, however, predicts large slip rates along the SWIM lineaments (up to 4 mm yr-1, depending on their length, Figs 8 and 9a), and little deformation in the northern Gulf of Cadiz and SW Iberia margin, in strong contrast with the distribution of strain computed from the instrumental seismicity (Figs 6 and S1e) and with the widespread evidences of recent re-activation of thrust faults to the north of the Horseshoe Abyssal Plain (e.g. Masson et al. 1994; Hayward et al. 1999; Terrinha et al. 2003; Zitellini et al. 2004; Terrinha et al. 2009). Moreover, the model is associated with a strong attenuation of the velocity field between northern Morocco and Gibraltar, in disagreement with present day GPS measurements (Fig. S1h). In our view, therefore, the SWIM lineaments should not be regarded as mature lithospheric-scale features marking the present day plate boundary between Nubia and Eurasia.
On the other hand, a tectonic model dominated by NNE-SSW to ENE-WSW thrust faults, locally disrupted by strike-slip (or transfer) faults (Model-3), is in particularly good agreement with the amount of shortening and wrenching between NW Morocco and SW Iberia estimated from GPS data (Fig. 12), the distribution of the seismic strain rate (Figs 11 and S1e), and the associated faulting mechanisms. In its widest segment, the region of transpressive deformation is over 200 km, between the Gulf of Cadiz and the northern Tagus Abyssal Plain (Fig. 2; Sartori et al. 1994; Masson et al. 1994; Hayward et al. 1999; Terrinha et al. 2003; Zitellini et al. 2004; Cunha et al. 2010), and extends over thinned continental, transitional and oceanic type crust of Jurassic-Early Cretaceous age (Gutscher et al. 2009, González et al. 1996, 2001; Rovere et al. 2004). Off SW Portugal, the model predicted fault slip rates are also consistent with the existing, though scarce, geological constraints, where the seismostratigraphy indicates a minimum thrusting rate of 0.25-0.5 mm yr-1 in the Marquês de Pombal Fault between the Tortonian (possibly base of Pliocene) and Recent (Zitellini et al. 2001; Roque 2008)
In the Gulf of Cadiz and SW Iberia margin, the convergence is accommodated by several structures, and maximum predicted long-term average faults slip rates vary between 1-2 mm yr-1; that is, <50 per cent the average plate relative movement. This region is known to have been the source of destructive earthquakes and tsunamis in the historical past, like the 1755 November 1 ‘Great Lisbon Earthquake’ with an estimated magnitude (Mw) between 8.5 and 8.7 (Johnston 1996; Solares & Arroyo 2004), and several other high-energy events since 7.0 kyr BP (c.f. Baptista & Miranda 2009 and references therein).
The low-predicted fault slip rates are thus an important result when assessing the seismic hazard of Portugal, SW Spain and Morocco, suggesting very long return periods for high-magnitude earthquakes on individual structures (exceeding a few thousand years). However, the existence of several mapped faults in the Gulf of Cadiz-SW Iberia margin that can generate Mw> 8 earthquakes (i.e. >60 km in length, approximately, Ribeiro et al. 2006; Stich , 2007; see tectonic map of Fig. 2 for location of the structures) or compound, multiple rupture hypothesis involving several fault systems (Zitellini et al. 2001; Terrinha et al. 2003; Ribeiro et al. 2006; Stich et al. 2007; Cunha et al. 2010; Terrinha et al. 2009), preclude a straightforward time dependent seismic and tsunami hazard assessment.
As discussed in Section 3.3, the crust and lithospheric structure used in the neotectonic models throughout this study, under the general assumptions of local isostasy and a steady state thermal regime, are broadly consistent with the observations and previous crust/lithosphere models in the region. The greater uncertainties are arguably in the eastern Gulf of Cadiz-Gharb basin-Rif area and under the Atlas mountains, where models such as those put forward by Zeyen et al. (2005) and Fullea et al. (2007, 2010) suggest a much thicker and thinner lithosphere, respectively.
We investigated the impact that such distinctive features may have in the modelling predictions by building an alternative lithospheric model, where we locally modified the heat flow input data (Fig. S3, Supporting Information). The model presented in Fig. S3(b) reproduces approximately the lithosphere in Fullea et al. (2007), with 140-160 km lithosphere thicknesses between the eastern Gulf of Cadiz and the Betic chain, extending into northern Morocco and a <80-km-thick lithosphere in Central Atlas. Despite the large differences in relation to the reference model (Fig. S3a), of up to 50 km in the Atlas mountains and 60-70 km in the Rif region, the predicted slip rates (Figs S3c and d), strain distribution and velocity field (Figs S3e and f) are very similar between the two models. The most notable effects are an increase in the amount of shortening in the Atlas mountains combined with a decrease in the thrusting rates along the Rif front. In terms of the scoring, the ‘Fullea-like’ lithosphere model is associated with a marginal (<5 per cent) decline in the correlation with the seismic strain rate, and a comparable improvement in the fit to the GPS data.
7 Conclusions
We summarize here the main results from a neotectonic modelling study of the western sector of the Nubia-Eurasia Plate boundary, between the Gloria Fault (off SW Iberia) and north Algeria (Fig. 1). In relation to previous neotectonic models in the region, we used a better constrained tectonic map for the SW Iberia and Gulf of Cadiz, and we utilize the recently made available GPS measurements to evaluate the modelling results, together with the instrumental seismicity and the stress data. In the SW Iberia margin and Gulf of Cadiz, three distinct tectonic models have been tested, based on the mapped structural elements, the seismicity of the area and proposed geodynamic models for region.
We show that the seismic strain rate, the stress field and the GPS velocities in Morocco and the Iberia Peninsula can be broadly explained assuming a relatively simple, two-plate tectonic framework, where Nubia and Eurasia converge NW-SE to WNW-ESE at a rate of 4.5-6 mm yr-1, as inferred from recent, geodetically constrained kinematic plate models. In relation to previous neotectonic models, covering both larger and more restricted areas of this sector of the Nubia-Eurasia Plate boundary, but which used the NUVEL-1A global plate model and the Argus 3-Plate model (Argus et al. 1989), an improvement of up to 40 per cent was achieved in terms of the correlation with the seismic strain rate.
In the Gulf of Cadiz and SW Iberia margin, the modelling results strongly support a tectonic model where the compressional deformation is accommodated along NNE-SSW to NE-SSW and ENE-WSW thrust faults and WNW-ESE right-lateral strike-slip faults. In our preferred model, ~1.8 mm yr-1 of shortening is predicted between northern Morocco and stable Nubia, ~1.4 mm yr-1 between NW Morocco and southernmost Portugal and ~0.4 mm yr-1 between southernmost Portugal and the latitude of Lisbon, associated with ~2.1 mm yr-1 of total right-lateral motion. These results are in excellent agreement with the GPS data. Moreover, the deviations between the model predicted and the GPS derived velocities around the Straits of Gibraltar and NW Morocco coast are <1 mm yr-1, which argues against the possibility of a second geodynamic engine to explain the data, such as an active subduction roll-back as proposed by Gutscher et al. (2002).
The proposed tectonic model is also consistent with the seismicity of the region, dominated by shallow to intermediate depth events with strike-slip and thrust focal mechanisms solutions, and the evidences of recent activity on numerous compressional and transpressional structures along the SW Iberia margin. The maximum predicted, long-term average fault slip rates vary between 1.5-2 mm yr-1 in the northern Gulf of Cadiz and 0.5-1 mm yr-1 in the NNE-SSW to NE-SW trending thrusts off SW Portugal, consistent with the existent geological constraints. The small slip rate values, only a fraction of the plate convergence velocity (4.5-6 mm yr-1), point to very long return periods for big earthquakes (Mw > 8) and tsunamis on individual structures.
On the other hand, a model which is dominated by long (up to 600 km), WNW-ESE trending lineaments across the Gulf of Cadiz and the Horseshoe abyssal plain (the SWIM lineaments) predicts little or no deformation in the faults located to the north, and a strong attenuation of the velocity field between the northern Morocco and the Straits of Gibraltar, in clear disagreement with the observations. The modelling results thus support a diffuse Nubia-Eurasia Plate boundary along SW Iberia and the Gulf of Cadiz, meaning that the emergent plate boundary along the recently mapped SWIM lineaments, as put forward by Zitellini et al. (2009), is still not mature and might have only a limited lithospheric expression.
Acknowledgments
The authors would like to gratefully acknowledge the funding from Fundação para a Ciência e Tecnologia (FCT, Portugal) through the research projects SWITNAME (PDCT/CGEGIN59244/2004), TOPOMED (TOPOEUROPE/0001/2007) and KINEMA (PTDC/CTEGIN/64101/2006), and the research fellowship SFRH/BPD/34798/2007. The authors would also like to acknowledge P. Bird and M. Sambridge for making the SHELLS code freely available, and are grateful to I. Jiménez-Munt for providing the initial heat-flow grid and help with the scoring of the seismic data. Constructive reviews by M. Miranda, Rob Govers and two anonymous reviewers helped improve an early version of this manuscript considerably.
References
