- Email: [email protected]
Crustal structure of the Betic–Rif system, western Mediterranean, from local earthquake tomography
Crustal structure of the Betic-Rif system, western Mediterranean, from local earthquake tomography Lahcen El Moudnib, Antonio Villase˜nor, Mimoun Harnafi, Josep Gallart, Antonio Pazos, Inmaculada Serrano, Diego C´ordoba, Javier A. Pulgar, Pedro Ibarra, Mohammed M. Himmi, Mimoun Chourak PII: DOI: Reference:
S0040-1951(15)00006-2 doi: 10.1016/j.tecto.2014.12.015 TECTO 126507
To appear in:
Tectonophysics
Received date: Revised date: Accepted date:
2 June 2014 14 November 2014 28 December 2014
Please cite this article as: El Moudnib, Lahcen, Villase˜ nor, Antonio, Harnafi, Mimoun, Gallart, Josep, Pazos, Antonio, Serrano, Inmaculada, C´ ordoba, Diego, Pulgar, Javier A., Ibarra, Pedro, Himmi, Mohammed M., Chourak, Mimoun, Crustal structure of the BeticRif system, western Mediterranean, from local earthquake tomography, Tectonophysics (2015), doi: 10.1016/j.tecto.2014.12.015
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT Crustal structure of the Betic-Rif system, western
SC RI
tomography
PT
Mediterranean, from local earthquake
Lahcen El Moudnib1,2 , Antonio Villaseñor3,* , Mimoun Harnafi1 , Josep
NU
Gallart3, Antonio Pazos4, Inmaculada Serrano5, Diego Córdoba6, Javier A.
MA
Pulgar7, Pedro Ibarra8, Mohammed M. Himmi2, and Mimoun Chourak9.
Institut Scientifique, Université Mohammed V-Agdal, Rabat, Morocco.
2
Faculté des Sciences, Université Mohammed V-Agdal, Rabat, Morocco.
3
Institute of Earth Sciences Jaume Almera, ICTJA-CSIC, Barcelona, Spain.
4
Real Instituto y Observatorio de la Armada, San Fernando, Cádiz, Spain.
5
Instituto Andaluz de Geofísica, Universidad de Granada, Granada, Spain.
6
IGEO and Universidad Complutense de Madrid, Madrid, Spain.
7
Departamento de Geodinámica Interna, Universidad de Oviedo, Oviedo, Spain.
8
Instituto Geológico y Minero de España, Madrid, Spain.
9
École Nationale des Sciences Appliquées, Université Mohammed Premier, Oujda,
AC
CE
PT
ED
1
Morocco.
*Corresponding author: [email protected]
Abstract We have determined the three-dimensional P-wave velocity structure of the crust and uppermost mantle beneath the Iberia-Africa collision zone using local earthquake
1
ACCEPTED MANUSCRIPT tomography. We have inverted arrival times of first-arriving P phases listed in the bulletins of the Instituto Geográfico Nacional and phases picked by us on permanent
PT
stations of other regional networks and on temporary broadband stations deployed in the frame of the TOPO-IBERIA project. In total we have used 38,624 P-wave arrival times
SC RI
from 2,362 local events recorded at 120 seismic stations. The most remarkable result is the imaging of a large high velocity body following approximately the western Alboran coastline, with a horizontal dimension of at least 200 km and extending in depth from
NU
the surface down to 24 km. This body, not imaged previously with this extent using
MA
seismic tomography, coincides with surface exposures of subcontinental mantle (peridotites) in Iberia and Africa and with a belt of positive gravity anomalies. We have also found a marked contrast in the seismic velocities of the middle and lower crust of
ED
the Alboran basin, coinciding with the location of the Trans-Alboran Shear Zone. We
PT
attribute this contrast to widespread magmatic intrusions in the eastern part of the basin, resulting in higher P-wave velocities than in the west. This contrast would also explain
CE
the different orientation of the Trans-Alboran Shear zone with respect to the surface
AC
features and faults in the Alboran basin. We also image thick crust beneath the Betics and Rif, accompanied by downgoing lithosphere of the Iberian foreland and Gulf of Cadiz beneath the Betic-Rif-Alboran system.
Keywords: tomography, Moho, Iberia, Morocco, Alboran
Highlights:
A high-velocity anomaly coincides with surface exposures of subcontinental mantle.
Different seismic structure between the western and eastern Alboran basin.
This velocity contrast might explain development of Trans-Alboran Shear Zone.
2
ACCEPTED MANUSCRIPT
The Iberian foreland basement deepens to the SE in the western Betics.
PT
1. Introduction The Iberia-Africa collision zone in the western Mediterranean is characterized by a
SC RI
complex interplate deformation zone. A tightly curved orogenic belt, the Betic-Rif system, is located between the undeformed Iberian and Morocco forelands (Figure 1). This section of the Eurasia-Africa plate boundary is characterized by a broad area of
NU
diffuse, moderate seismicity in contrast with the adjacent segments to the east (Algeria)
MA
and to the west (Gulf of Cadiz) where large earthquakes (M > 7) are frequent. The Betic-Rif system is traditionally divided into internal and external zones (e.g. Platt et al.,
ED
2013 and references therein). The external zones consist primarily of folded Mesozoic and Cenozoic sedimentary rocks from the former Iberian and African paleomargins,
PT
while the internal zones are composed of a stack of pre-Mesozoic and Mesozoic metamorphic rocks. Within the internal zones there are several exposures of lithospheric
CE
mantle with varying sizes: the Ronda peridotite in southern Spain is the largest outcrop
AC
of subcontinental mantle in the world with an extent up to 300 km2 (Obata, 1980); the Beni Bousera peridotite, with 70 km2 is the equivalent outcrop in the Rifean range; and finally some minor outcrops are also found in the Sarchal beach in Ceuta (SánchezGómez et al., 1995). The Alboran basin, located within the internal part of the Betic-Rif orogen, is an extensional basin formed during the Neogene (Comas et al., 1999). It is floored in the west by continental basement that largely corresponds to rocks belonging to the internal complexes of the Betics-Rif and by magmatic arc crust in the east (Booth-Rea et al., 2007). Although the geology and shallow structure of the region have been intensively studied (see recent reviews by Gutscher et al., 2012; Platt et al., 2013) the debate over the
3
ACCEPTED MANUSCRIPT kinematics and dynamics of its formation is still intense. Alternative evolution models based on kinematic reconstructions have been proposed for the region by Faccenna et
PT
al. (2004), Vergés and Fernández (2012) and van Hinsbergen et al. (2014) to mention only a few. Chertova et al. (2014) have tested these models using 4D thermomechanical
SC RI
numerical modeling and have been able to reproduce the slab morphology in the mantle observed using teleseismic tomography.
With respect to the deep structure of the crust and upper mantle, the first order features
NU
have been know for some time, for example the presence of a large high velocity
MA
anomaly in the upper mantle beneath the Alboran basin that probably corresponds to a subducted slab (Blanco and Spakman, 1993). Similarly, a significant number of studies have used local earthquake tomography with arrival times recorded by local networks to
ED
image the crustal structure of the region (Dañobeitia et al., 1998; Gurria and Mezcua,
PT
2000; Morales et al., 1999; Serrano et al., 2002; Serrano et al., 2003; Timoulali et al., 2014) . However, many of them are based on arrival times picked on analog recordings,
CE
and their resolution in the Alboran basin and northern Morocco is limited because of the
AC
sparse station coverage. In this study we make use of new high-quality data provided by the digitally-upgraded local seismic networks and also by a dense deployment of temporary broadband stations in southern Spain and northern Morocco. With this improved dataset we aim at obtaining more detailed images of the crust and uppermost mantle of the Iberia-Africa collision zone that can provide additional constraints on its past and present geodynamic evolution.
2. Data The data used in this study are first arrival times of P phases of local earthquakes, recorded at seismic stations. Because of the characteristics of the methodology, only
4
ACCEPTED MANUSCRIPT earthquakes and stations that are located inside the study region (model volume) are used. Our main source for arrival times is the bulletins of the Instituto Geográfico
PT
Nacional (IGN) from 2000-2009. The IGN operates the Spanish National Seismic Network (RSN, Red Sísmica Nacional), which by the year 2000 had already completed
SC RI
its upgrade to broadband sensors, digital recordings, and satellite communications. Therefore the phase arrival time data after 2000 are presumably of higher quality than earlier analog data. The IGN bulletins contain not only arrival times from stations of the
NU
Spanish RSN, but also include arrival times from neighboring networks such as the
MA
national Portuguese and Moroccan monitoring networks, operated by the Instituto Português do Mar e da Atmosfera (IPMA) and Centre National pour la Recherche Scientifique et Technique (CNRST) respectively. Additionally they include arrival time
ED
data from other research networks in the region, for example those operated by the
PT
Spanish Royal Navy Observatory (ROA). The station coverage provided by these networks is well suited for locating local earthquakes in the southern Iberian Peninsula,
CE
because of the greater station density in the most seismically active onshore areas.
AC
However this situation is less optimal for seismic tomography studies, for which a homogeneous station coverage would be preferable. The situation is still worse for northern Morocco because of the larger station spacing there and the lower quality of the arrival time data (for the time period considered here the national network was still analog). To ameliorate this situation the TOPO-IBERIA project, launched in 2007, included the deployment of the IberArray temporary network (Díaz et al., 2009; Díaz et al., 2010) with the objective of achieving a homogeneous station coverage in Spain with an interstation distance of 60 km or smaller. Simultaneously 18 temporary broadband stations with the same instrumentation as IberArray were also installed in northern
5
ACCEPTED MANUSCRIPT Morocco, in order to reduce the low station coverage provided by permanent networks. Subsequently other temporary networks (e.g. PICASSO, WILAS) have been installed in
PT
the Iberian Peninsula and Morocco both expanding and densifying the coverage of
SC RI
IberArray, but their data have not been incorporated in this study.
We have complemented the arrival time database of the IGN bulletins from 2000 to 2009 with phase picks made by us on digital recordings from the IberArray stations in
NU
Spain and Morocco for selected earthquakes in the time period 2007-2009. We have
MA
also incorporated picks made by us on recordings from other permanent networks that are not routinely included in the IGN bulletins, such as the Andalusian Seismic Network operated by the Instituto Andaluz de Geofísica (IAG). In total we have used data from
ED
120 stations in the region, 55 of them new broadband instruments previously
PT
unavailable (35 in the Iberian Peninsula and 20 in northern Africa). Figure 2 illustrates the improvement in coverage and quality provided by the new dataset. The use of a
CE
seismic array with a dense, homogeneous spacing not only improves the data coverage
AC
but also allows to identify errors, for example stations with residuals significantly different in size and/or sign than those of neighboring stations. In most cases these differences can be traced to bad pickings due to low signal-to-noise ratio or timing errors (stations that temporary lost the GPS time signal). The use of a seismic array also allows to evaluate qualitatively the strength and patterns present in the data. For example the distribution of residuals for the NE Morocco earthquake shown in Figure 2a clearly indicates a correlation with the different structural units. The Iberian and Moroccan forelands (Iberian Massif and Moroccan Meseta) and the Betics internal zones are clearly associated with negative (fast) residuals, while the external Betics and Rif exhibit large positive (slow) residuals. This general pattern is also observed for the
6
ACCEPTED MANUSCRIPT Gulf of Cadiz event shown in Figure 2b, although the negative residuals in the internal Betics and Portuguese Iberian Massif are now smaller in size or even positive. Stations
PT
located next to the Spanish Alboran shore are almost always consistently negative, indicating the presence of thinner than average crust and/or the presence of a large high-
SC RI
velocity body.
From this combined dataset we selected well recorded earthquakes for the tomographic
NU
inversion. First we choose a rectangular region of 880 720 km2 centered in the
MA
Alboran basin with a good distribution of local earthquakes and seismic stations (map in Figure 3 corresponds to the target region). Seismicity in this region is characterized by
ED
large variability in focal depth, including shallow, intermediate-depth and even a cluster of deep-focus earthquakes at approximately 650 km depth. Because of the decrease in
PT
seismicity with depth and in order to keep the model size manageable, we have considered only earthquakes with focal depth h < 100 km. We have then extracted
CE
earthquakes inside this volume with P phases recorded by 10 or more stations at
AC
distances less than 300 km, and with an azimuthal gap smaller than 200. Koulakov (2009) has shown that the use of events with azimuthal gaps greater than the commonly used value of 180 can contribute to improving the images obtained from local earthquake tomography. Here we use a value only slightly larger than 180 in order to try to improve the ray coverage in the marine parts of the model of interest, particularly in the Gulf of Cadiz. We have also required that selected earthquakes include 2 or more S readings because of their positive influence in improving the determination of focal depths (Gomberg et al., 1990). Finally, to avoid large outliers and misidentifications, we have eliminated earthquakes whose hypocenter solution reported by IGN had root mean square (RMS) residuals larger than 3.5 s. These selection criteria resulted in a dataset of
7
ACCEPTED MANUSCRIPT 38,624 P-wave arrival times from 2,362 local events recorded by 120 seismic stations, shown in Figure 3. Selected seismicity is rather well distributed over the entire model
PT
region, but three prominent seismic zones can be distinguished. The Trans-Alboran Shear Zone (TASZ) is a linear feature that extends approximately from Al-Hoceima in
SC RI
Morocco to southeastern Spain. A dense cluster of seismicity in the western Betics near the town of Morón is due to two moderate events (Mw ~ 4.5) in 2007-2008 and their aftershock sequences. While the seismicity in the region is mostly shallow (crustal), a
NU
remarkable zone of subcrustal, intermediate-depth seismicity occurs in the western
MA
Alboran basin, extending approximately N-S from Malaga to the west of Al-Hoceima. These deep earthquakes contribute to the dataset with ray-paths that improve the
ED
illumination of the uppermost mantle.
PT
3. Method
The tomographic inversion method used in this study is described in detail by Benz et
CE
al. (1996) and will only be summarized here. It has been successfully applied to
AC
different spatial scales, including volcanoes (Benz et al., 1996; Okubo et al., 1997; Villasenor et al., 1998) and to larger regions such as the one considered here (e.g. Powell et al., 2014). In this technique, arrival times of first-arriving P waves are inverted to simultaneously determine the three-dimensional P-wave velocity structure and earthquake relocations. This non-linear problem is linearized and solved iteratively. In each iteration, updated velocity model and earthquake locations are obtained and used as starting values for the next iteration. The iteration procedure is stopped after a preset number of iterations (in our case 10 iterations) or when the data misfit is not reduced significantly. Theoretical travel times are calculated using the finite-difference technique of (Podvin and Lecomte, 1991), which is well suited for computing travel
8
ACCEPTED MANUSCRIPT times in media with large lateral velocity perturbations. Velocity anomalies with respect to a reference model are solved using the iterative LSQR method of Paige and Saunders
PT
(1982) with smoothing constraints applied to control the degree of model roughness. We have tested a wide range of values of the smoothing parameter k (see Equation 7 in
SC RI
Benz et al. (1996) for its definition) in order to find those that produce good image fidelity without introducing large velocity alternations in adjacent cells and at the same time produce a significant reduction in data misfit. We have tested values of k ranging
NU
from 300 (smooth model) to 10 (very rough). Results show that decreasing the
MA
smoothing parameter down to a value of 75 continuously improves the data fit. Smaller values result in very rough models, and oscillations in the data fit with each iteration. Therefore we have adopted a value of k = 75, because it is the one that produces the best
ED
variance reduction without degrading the appearance of the model. A more detailed
PT
description of these tests is shown in the Supplementary Material, together with examples of under- and over-damped models.
CE
The reference/starting model considered here is the 1-D layered model used by the IGN
AC
for routine earthquake location in the Iberian Peninsula (Table 1). This model is representative of the velocity structure of the Iberian Massif (Paleozoic units in Figure 1) and has a crustal thickness of 31 km. We have interpolated the constant-velocity layers and discontinuities of the IGN 1-D model to linear gradients with smooth variations, in order to prevent that rays travel exclusively along the layer discontinuities and do not sample other parts of the model. However, if the arrival time data contains enough paths refracted through subsurface discontinuities (e.g. Pn phases) these discontinuities will naturally develop in the inverted tomographic model. The model volume, as described in the previous section, has horizontal dimensions of 880 720 km2 and, in order to include all potential ray-paths between earthquakes and
9
ACCEPTED MANUSCRIPT stations, it extends from 4 km above sea level to 120 km depth. The velocity model is parameterized as a regular grid of constant velocity cells of 20 20 4 km3 (in the x, y
PT
and z directions respectively), resulting in a total of 49,104 model parameters. A finer grid of 1 1 1 km3 was used for the finite-difference traveltime calculations. The size
SC RI
of the model and traveltime cells was determined empirically through a series of tests that addressed the effects of cell size and smoothing constraints (Benz et al., 1996).
NU
4. Results
MA
Inversion of the P-wave arrival time data after 10 iterations resulted in a reduction of 28% of the residual RMS from an starting value of 0.83 s. Because of the large number
ED
of unknowns, we cannot compute the formal model resolution and covariance, and therefore the resolution analysis of our model is done empirically using the ray-path
PT
coverage and synthetic reconstruction checks using the source-receiver geometry of our dataset. We use spike tests, that are similar to traditional checkerboard tests, except that
CE
the alternating positive and negative anomalies are not adjacent, but separated by a
AC
buffer with 0% velocity anomaly. These results show good recovery of the synthetic anomalies without significant horizontal smearing in the central part of the model, which is well sampled by through-going rays (Figure 4a,b). Resolution decreases below 40 km where most rays have vertical paths coming from intermediate-depth earthquakes in the western Alboran region. Hit count maps and results from spike tests for two layers of the model are shown in Figure 4 and with more detail in the Supplementary Material. As a result of the simultaneous inversion for velocity structure and earthquake location, new locations have been obtained. However, changes in location are not very large. The average difference in epicenter between the original IGN locations and the final
10
ACCEPTED MANUSCRIPT relocations is 4.5 km, and the average difference in depth is -1 km (the minus sign indicates that the final locations are shallower on average than the initial ones).
PT
In order to evaluate the influence of the added arrival times measured on IberArray stations in southern Spain and northern Morocco we have performed tomographic
SC RI
inversions with and without including this dataset. The result of this comparison is shown in Figure 5. The top panel shows a horizontal slice of the model (8-12 km depth) obtained using data listed in the IGN bulletins alone, while the bottom panel shows the
NU
same horizontal slice obtained using the combined dataset. An obvious improvement is
MA
the extension to the south of the region that is well covered by ray paths to include most of the Rif system in northern Morocco. There are, however some substantial differences even in regions that are reasonably well sampled by the IGN dataset. The most
ED
remarkable of these differences is in the shape and extension of the high velocity
PT
anomaly associated with the outcrop of the Ronda peridotites. While in the results using only the IGN bulletins this anomaly appears as a large blob, in the results from the
CE
combined dataset it has a curved, narrow, elongated shape, coincident with the Iberian
AC
coastline of the Alboran basin, which is in good agreement with the distribution of positive residual gravity anomalies in the region (Torne et al., 2000). Another relevant difference is the sharp contrast in velocity observed in the results of the combined dataset across the TASZ, with very low velocities in the west and faster velocities in the east, which is absent in the results based on the IGN bulletins alone. Figure 6 shows other horizontal layers of the model from the shallow crust to uppermost mantle. A common feature to all layers shown is the large variability in the values of the velocity perturbation with respect to the reference model, which often exceed 10%. This is particularly true for the negative (slow) velocity anomalies, which extend through large regions of the model.
11
ACCEPTED MANUSCRIPT Sedimentary basins are generally characterized by large low velocity anomalies. However, in our model the most important sedimentary basins (Guadalquivir, Gharb,
PT
and West Alboran basins) are poorly imaged because few stations and earthquakes are located on them, resulting in a lack of sampling rays. This is not the case for the
SC RI
intramontane basins in the Betic range, that are characterized by moderate seismicity and well monitored by seismic networks. The Granada-Guadix-Baza basins can be clearly identified as low velocities in the shallow layers of the model (Figure 6a).
NU
Similarly, the land portion of the Gulf of Cadiz accretionary complex is characterized
MA
by slow velocities, which are remarkably low in southernmost Spain, coinciding with the location of the Campo de Gibraltar flyschs. These low velocities not only occur at shallow depths, where they can be associated with the flysch units, but throughout the
ED
entire crust. Another conspicuous low velocity anomaly is located in the central Alboran
PT
basin and follows the trend of the TASZ. Although small sedimentary basins occur in this region, this anomaly is too large to be explained exclusively by sediments.
CE
The curved high velocity anomaly associated with the surface expression of the Ronda
AC
peridotite can be clearly followed from shallow depths (e.g. Figure 5b) to 12-16 km depth (Figure 6b) and down to at least 24 km (not shown here). This is also the case for the high velocity anomaly located east of the TASZ that is clearly imaged from 8 to 24 km depth. Figure 6c shows the layer of the model between 28 and 32 km depth, which is located just above the Moho of the reference model. A conspicuous feature in this layer is a high velocity zone that coincides with the extent of the Alboran basin. This is most likely caused by crustal thickness in this region being smaller than that of the reference model (31 km). Therefore the mantle at this depth would have higher velocities than the surrounding crustal blocks. The Ronda peridotite and eastern TASZ high velocity
12
ACCEPTED MANUSCRIPT anomalies described previously seem to be rooted in this high velocity anomaly,
PT
indicating that they extend throughout the crust down to the top of the mantle.
Figure 6d shows the layer between 36 and 40 km depth, within the mantle of the
SC RI
reference model. Large low velocity anomalies in this layer could indicate crustal velocities and therefore crustal thicknesses larger than the reference model. This is in fact observed all along the Betic-Rif system and would indicate the presence of a crustal
NU
root beneath those ranges.
MA
5. Discussion
ED
Here we will discuss the most relevant results obtained in this study and their
PT
implication for the geodynamic evolution of the region and present-day tectonics.
Anomalies associated with peridotite outcrops
CE
The most remarkable new finding in this study is the curved, narrow, high-velocity
AC
anomaly that follows the northwestern shore of the Alboran basin and coincides with the exposure of subcontinental mantle in the Ronda peridotite. Other studies based on local earthquake (Timoulali et al., 2014) and surface wave data (Chourak et al., 2005) have found high seismic velocities coinciding with the Ronda peridotites, but not with the geometry and geographical extent obtained here. This can be attributed to the lack of ray paths to stations in Morocco, and underscores the importance of the new dataset provided by the IberArray temporary stations. In our model this anomaly can be traced continuously for at least 200 km from approximately 3.5W in the east up to Ceuta in the southern shore of the Strait of Gibraltar. (e.g. Figure 6b), suggesting that the Ronda and Ceuta peridotite outcrops correspond to the same continuous body at depth.
13
ACCEPTED MANUSCRIPT Although the high velocity anomaly imaged in our model continues south of Ceuta it does not reach unequivocally the location of the Beni Bousera peridotites (Figure 6b).
PT
However other geophysical evidence such as gravity anomalies (Torne et al., 2000) and ambient noise tomography (Silveira et al., 2013; Villaseñor et al., 2011) suggest that
SC RI
this anomaly extends further south and therefore encompasses the three exposures of subcontinental mantle in the region. The inability of our model to image the continuation of the high velocity into the Beni Bousera peridotite is probably due to the
NU
lower station density in northern Morocco. As seen in Figure 1, few stations (none of
MA
them broadband) are located next to the Moroccan western Alboran shore nor in the Rif internal zones, preventing the detailed imaging of the structure in that region. Incorporation of other existing stations (such as those from the PICASSO experiment)
ED
would be necessary to improve the resolving power of the traveltime dataset in the
PT
internal Rif.
In the same way this high velocity anomaly has a large lateral extent, it also exhibits a
CE
large vertical continuity, extending from the surface to at least 24 km depth. This can be
AC
clearly observed in the vertical cross section of the model shown in Figure 7a. This figure also shows that the anomaly is not vertical but dips towards the center of the Alboran basin.
The Trans-Alboran Shear Zone (TASZ) Another new remarkable feature of the model is the strong change in P-wave velocity across the TASZ, with slow velocities in the west and fast in the east. This contrast, however is not observed at all depths, but most clearly between 8 and 24 km. The location of this boundary can be seen in the vertical cross section shown in Figure 7b. The orientation of the TASZ (N35E, Stich et al., 2006) is somewhat unexpected because
14
ACCEPTED MANUSCRIPT it is nearly perpendicular to the direction of convergence between Eurasia and Africa at this location, and is also different to the SW-NE orientation of the main sinistral faults
PT
in the Rif (Jebha and Nekor) and bathymetric features in the Alboran basin (e.g. the Alboran ridge). However, if the high velocity anomalies east of the TASZ correspond to
SC RI
a mechanically stronger block of the Alboran basin crust, it could act as an indenter causing the extrusion to the SW of a weaker western Alboran crust, thus explaining the left-lateral strike-slip mechanisms along the TASZ (see for example Figure 7 in Stich et
NU
al., 2003). Since the high velocity body is only observed below 8 km depth, it is
MA
compatible with structures in the upper crust having different orientations, for example inherited from the extension episode that generated the Alboran basin. The origin of this strong, high-velocity block, which coincides approximately with the eastern Alboran
ED
volcanic belt drawn by Duggen et al. (2003), could be magmatic intrusions in the
PT
preexisting thinned continental crust. In fact, magmatic arc crust is observed in the eastern Alboran basin, and arc-magmatism intrusions are more common in the eastern
CE
Alboran basin than in the west (Booth-Rea et al., 2007). The boundaries of crustal types
AC
determined by Booth-Rea et al. (2007) do not coincide exactly with the boundary of the block imaged here east of the TASZ but since those are based on shallow seismic studies, it is not incompatible with a different geometry at depth. The high velocity block seems to extend inland in Morocco (Figures 5b and 6b), in general agreement with previous tomographic studies (Serrano et al., 2003). In our model the boundary between high and low velocities seems to be located just east of Al-Hoceima, coinciding with an exposure of arc volcanic rocks. The inland extension of the western boundary of the high velocity block is also coincident with a region of a pronounced change in crustal thickness, with values of 25 km in the east to 40 km and deeper to the west (Mancilla et al., 2012).
15
ACCEPTED MANUSCRIPT
The onshore Gulf of Cadiz accretionary complex
PT
Very low velocities are observed at shallow depths in the westernmost external Betics. This anomaly coincides with the largest outcrop of flysch units in the region (e.g.
SC RI
Crespo-Blanc et al., 2012). Interestingly at all crustal depths of our model, this region is the one than exhibits the lowest seismic velocities. This observation has also been found independently by Silveira et al. (2013) using seismic ambient noise. In their study they
NU
obtain very low velocities (up to 40% slower than the average) for the propagation of
MA
Rayleigh waves travelling through the crust in the westernmost Betics and Rif. This indicates that P-wave velocities in the western Rif should be similarly slow to those in the western Betics if a comparable station and ray-path coverage were available.
ED
These anomalously low P- and surface-wave velocity values suggest that the flysch
PT
trough units and/or sediments of the Iberian-Morocco paleomargins have been underthrusted to great depths. This has been suggested for the Rif by Chalouan et al.
CE
(2001) and is also suggested by the cross sections in Figure 7 for the western Betics.
AC
Figure 7a clearly shows low-velocity rocks overridden by the high velocity anomaly associated to the Ronda peridotite. In map view there is also a good correlation between the location of this high velocity anomaly and the exposures of the flysch trough units (Figure 1) suggesting that the high velocity anomaly might have acted as a backstop for the underthrusting of the flyschs and other units.
The Iberian foreland basement Incorporation in the IGN bulletins of stations located in southern Portugal provide some constraints on the seismic structure of the southwestern Iberian Massif. Since the reference model used in this study (1-D model used by IGN for earthquake location) is
16
ACCEPTED MANUSCRIPT representative of the structure in the Iberian Massif, velocity anomalies in that region should be close to zero. This is indeed the case for the most of the crustal layers
PT
(Figures 6a, 5b and 6b). However, when we approach to the crust-mantle discontinuity some differences begin to emerge. In Figure 6c, corresponding to the layer of the model
SC RI
representing the reference Moho, we observe a substantial difference between the western and eastern Iberian Massif. Beneath Portugal the velocity anomalies continue to be close to zero or slightly negative. This behavior is already clearly present in the
NU
arrival time data (see the residual pattern in the Gulf of Cadiz earthquake shown in
MA
Figure 2b). Conversely, velocities in the eastern Iberian Massif are larger than average and extend farther south into the deformation front of the external Betics. It is not clear what is the reason for this east-west change in seismic velocities, because it does not
ED
coincide with any of the boundaries between the units that form the Iberian Massif.
PT
Crustal thickness throughout the Iberian Massif is remarkably constant, approximately 31 km (e.g. Díaz and Gallart, 2009), and this can be observed in the northern part of the
CE
cross section shown in Figure 7a. However, towards the SE there is a clear deepening of
AC
the contrast between crustal and mantle velocities. This deepening is also accompanied by a SE dipping zone of intermediate-depth earthquakes, that can be followed up to the surface. Similarly, although not as well constrained as for the Iberian Massif, the crustmantle discontinuity beneath the Gulf of Cadiz and the earthquake focal depths seem to dip towards the east (Figure 7b). This suggest that the continental lithosphere of Iberia (and also the lithosphere beneath the Gulf of Cadiz) is undergoing continental subduction (or alternatively delamination or roll-back) beneath the Alboran Sea. Furthermore, in this part of the western Betics, the Iberian lithosphere would still be attached to the slab that has been imaged in the mantle beneath the Alboran Sea (e.g. Garcia-Castellanos and Villaseñor, 2011; Spakman and Wortel, 2004). A similar
17
ACCEPTED MANUSCRIPT observation has been obtained using receiver functions (Mancilla et al., 2013). This process of continental subduction/delamination could also be present in the Rif, but
SC RI
Crustal thickness beneath the Betics and Rif
PT
cannot be constrained with our present dataset.
Regions in which the crust is thicker than the reference Moho would be characterized in Figure 6c by significant low velocity anomalies, which correspond exactly to the entire
NU
Betic-Rif system, including areas such as the Strait of Gibraltar that are not
MA
characterized by high topography. A very similar pattern in crustal thickness is also observed using receiver functions (Thurner et al., 2014; Mancilla et al., submitted), and ambient noise tomography (Villaseñor et al., 2011). Between 36-40 km depth (Figure
ED
6d) strong low velocities are still observed beneath the central Betics, indicating the
CE
6. Conclusions
PT
existence of a deeper root beneath the areas of highest topography.
AC
We have obtained a new P-wave velocity model of the Iberia-Africa collision zone, that includes the Betic-Rif system and the Alboran basin. A major finding of this study is a large high velocity body 200-km long and extending from 0 to approximately 24 km depth that is associated with the surface exposures of subcontinental mantle in the Ronda and Sarchal (Ceuta) peridotites, and possibly extending to the Beni Bousera peridotites in the Rif. Because of its size and location, this body, not previously imaged with this extent by similar studies, might hold important clues with respect to the evolution of the Alboran-Betic-Rif region. Likewise, because of its potential effect in teleseismic travel times and surface wave dispersion, it should be taken into account for seismic studies of the mantle beneath this region that do not consider crustal structure
18
ACCEPTED MANUSCRIPT and use instead crustal corrections. We have also imaged the different character of the eastern and western Alboran middle and lower crust, and its potential implications for
PT
the development of the TASZ. Most of these results have been possible because of the incorporation of arrival times
SC RI
measured on the temporary IberArray broadband stations in southern Spain and Morocco. Due to the lower station density in northern Morocco (although greatly improved with respect to previous studies) some features such as the continuation of the
NU
high velocity body into the Beni Bousera peridotite, the estimate of crustal thickness
MA
beneath the Rif, and the behavior of the African foreland basement towards the Alboran basin have not been completely resolved. Incorporation of data from other existing temporary deployments in the region (both marine and on land) might help to resolve
PT
ED
these issues.
CE
Acknowledgements
(IGN, Spain) for providing the bulletin and
AC
waveform data that are the base of this study. We also thank other networks operating in the region that have provided waveform data for this study: Instituto Andaluz de Geofísica (IAG, Spain), Real Instituto y Observatorio de la Armada (ROA, Spain), Instituto Português do Mar e da Atmosfera (IPMA, Portugal). Data from other individual stations have been obtained from the ORFEUS (The Netherlands), GEOFON (Germany) and IRIS (USA) data centers. This is a contribution of the Team ConsoliderIngenio 2010 TOPO-IBERIA (CSD2006-00041). Additional funding was provided by the SIBERIA (CGL2006-01171), RIFSIS (CGL2009-09727) and ALERTES-RIM (CGL2013-45724-C3-3-R) projects.
19
ACCEPTED MANUSCRIPT References Benz, H.M., Chouet, B.A., Dawson, P., Lahr, J.C., Page, R.A., Hole, J.A., 1996. Threedimensional P and S wave velocity structure of Redoubt Volcano, Alaska J.
PT
Geophys. Res. 101, 8111-8128.
Blanco, M.J., Spakman, W., 1993. The P-wave velocity structure of the mantle below
Tectonophysics 221, 13-34. -
-
SC RI
the Iberian Peninsula: evidence for subducted lithosphere below southern Spain.
, J.M., Grevemeyer, I., 2007. Crustal
types and Tertiary tectonic evolution of the Alborán sea, western Mediterranean.
NU
Geochem. Geophys. Geosyst. 8, Q10005.
Chalouan, A., Michard, A., Feinberg, H., Montigny, R., Saddiqi, O., 2001. The Rif
MA
mountain building (Morocco): a new tectonic scenario. Bulletin de la Société Géologique de France 172, 603-616.
Chertova, M.V., Spakman, W., Geenen, T., van den Berg, A.P., van Hinsbergen, D.J.J.,
ED
2014. Underpinning tectonic reconstructions of the western Mediterranean region with dynamic slab evolution from 3-D numerical modeling. Journal of Geophysical
PT
Research: Solid Earth 119, 5876-5902. Chourak, M., Corchete, V., Badal, J., Gómez, F., Serón, J., 2005. Shallow Seismic
CE
Velocity Structure of the Betic Cordillera (Southern Spain) from Modelling of Rayleigh Wave Dispersion. Surv. Geophys. 26, 481-504.
AC
Comas, M.C., Platt, J.P., Soto, J.I., Watts, A.B., 1999. The origin and tectonic history of the Alboran Basin: insights from LEG 161 results, in: Zahn, R., Comas, M.C., Klaus, A. (Eds.). Ocean Drilling Program, College Station, TX, USA, pp. 555-580. Crespo-Blanc, A., Balanyá, J.C., Expósito, I., Luján, M., Suades, E., 2012. Crescentlike large-scale structures in the external zones of the western Gibraltar Arc (BeticRif orogenic wedge). Journal of the Geological Society 169, 667-679. Dañobeitia, J.J., Sallarès, V., Gallart, J., 1998. Local earthquakes seismic tomography in the Betic Cordillera (southern Spain). Earth Planet. Sci. Lett. 160, 225-239. Díaz, J., Gallart, J., 2009. Crustal structure beneath the Iberian Peninsula and surrounding waters: A new compilation of deep seismic sounding results. Phys. Earth Planet. Inter. 173, 181-190. Díaz, J., Villaseñor, A., Gallart, J., Morales, J., Pazos, A., Córdoba, D., Pulgar, J.A., García-Lobón, J.L., Harnafi, M., Group, T.-I.S.W., 2009. The IBERARRAY
20
ACCEPTED MANUSCRIPT broadband seismic network: a new tool to investigate the deep structure beneath Iberia. Orfeus Newsl. 8, 1-6. Díaz, J., Villaseñor, A., Morales, J., Pazos, A., Córdoba, D., Pulgar, J.A., Garcia-Lobon,
PT
J.L., Harnafi, M., Carbonell, R., Gallart, J., 2010. Background Noise Characteristics at the IberArray Broadband Seismic Network. Bull. Seism. Soc. Am. 100, 618-628.
SC RI
Duggen, S., Hoernle, K., van den Bogaard, P., Rüpke, L.H., Phipps Morgan, J., 2003. Deep roots of the Messinian salinity crisis. Nature 422, 602-606. Faccenna, C., Piromallo, C., Crespo-Blanc, A., Jolivet, L., Rosseti, F., 2004. Lateral slab deformation and the origin of the western Mediterranean arcs. Tectonics 23,
NU
TC1012.
Garcia-Castellanos, D., Villaseñor, A., 2011. Messinian salinity crisis regulated by
MA
competing tectonics and erosion at the Gibraltar arc. Nature 480, 359-363. Gomberg, J.S., Shedlock, K.M., Roecker, S.W., 1990. The effect of S-wave arrival times on the accuracy of hypocenter estimation. Bull. Seism. Soc. Am. 80, 1605-
ED
1628.
Gurria, E., Mezcua, J., 2000. Seismic tomography of the crust and lithospheric mantle
PT
in the Betic Cordillera and Alboran Sea. Tectonophysics 329, 99-119. Gutscher, M.A., Dominguez, S., Westbrook, G.K., Le Roy, P., Rosas, F., Duarte, J.C.,
CE
Terrinha, P., Miranda, J.M., Graindorge, D., Gailler, A., Sallares, V., Bartolome, R., 2012. The Gibraltar subduction: A decade of new geophysical data. Tectonophysics 574-575, 72-91.
AC
Koulakov, I., 2009. Out-of-Network Events Can Be of Great Importance for Improving Results of Local Earthquake Tomography. Bull. Seism. Soc. Am. 99, 2556-2563. Mancilla, F., Stich, D., Berrocoso, M., Martin, R., Morales, J., Fernández-Ros, A., Paez, R., Pérez-Peña, A., 2013. Delamination in the Betic Range: Deep structure, seismicity, and GPS motion. Geology 41, 307-310. Mancilla, F.d.L., Stich, D., Morales, J., Julià, J., Díaz, J., Pazos, A., Córdoba, D., Pulgar, J.A., Ibarra, P., Harnafi, M., Gonzalez-Lodeiro, F., 2012. Crustal thickness variations in northern Morocco. J. Geophys. Res. 117, B02312. Morales, J., Serrano, I., Jabaloy, A., Galindo-Zaldívar, J., Zhao, D., Torcal, F., Vidal, F., González-Lodeiro, F., 1999. Active continental subduction beneath the Betic Cordillera and the Albor
. Geology 27, 735-738.
21
ACCEPTED MANUSCRIPT Obata, M., 1980. The Ronda Peridotite: Garnet-, Spinel-, and Plagioclase-Lherzolite Facies and the P-T Trajectories of a High-Temperature Mantle Intrusion. J. Petrol. 21, 533-572.
PT
Okubo, P.G., Benz, H.M., Chouet, B.A., 1997. Imaging the crustal magma sources beneath Mauna Loa and Kilauea volcanoes, Hawaii. Geology 25, 867-870.
SC RI
Paige, C.C., Saunders, M.A., 1982. LSQR: An algorithm for sparse linear equations and sparse least squares. ACM Trans. Math. Software 8, 43-71.
Platt, J.P., Behr, W.M., Johanesen, K., Williams, J.R., 2013. The Betic-Rif Arc and Its Orogenic Hinterland: A Review. Annu. Rev. Earth Planet. Sci. 41, 313-357.
NU
Podvin, P., Lecomte, I., 1991. Finite difference computation of traveltimes in very contrasted velocity models: a massively parallel approach and its associated tools.
MA
Geophys. J. Int. 105, 271-284.
Powell, C.A., Withers, M.M., Cox, R.T., Vlahovic, G., Arroucau, P., 2014. Crustal velocity structure associated with the eastern Tennessee seismic zone: Vp and Vs
ED
images based upon local earthquake tomography. Journal of Geophysical Research: Solid Earth 119, 464-489.
PT
Sánchez-Gómez, M., García Dueñas, V., Muñoz, M., Balanyá, J.C., 1995. Relación estructural entre los cuerpos peridotíticos situados al Norte y al Sur del Estrecho de
CE
Gibraltar. Geogaceta 17, 135-137. Serrano, I., Zhao, D., Morales, J., 2002. 3-D crustal structure of the extensional Granada Basin in the convergent boundary between the Eurasian and African plates.
AC
Tectonophysics 344, 61-79. Serrano, I., Zhao, D., Morales, J., Torcal, F., 2003. Seismic tomography from local crustal earthquakes beneath eastern Rif Mountains of Morocco. Tectonophysics 367, 187-201. Silveira, G., Dias, N.A., Villaseñor, A., 2013. Seismic imaging of the western Iberian crust using ambient noise: Boundaries and internal structure of the Iberian Massif. Tectonophysics 589, 186-194. Spakman, W., Wortel, M.J.R., 2004. A tomographic view on Western Mediterranean Geodynamics, in: Cavazza, W., Roure, F., Spakman, W., Stampfli, G.M., Ziegler, P.A. (Eds.), The TRANSMED Atlas, The Mediterranean Region from Crust to Mantle, pp. 31-52. Stich, D., Ammon, C.J., Morales, J., 2003. Moment tensor solutions for small and moderate earthquakes in the Ibero-Maghreb region. J. Geophys. Res. 108, 2148. 22
ACCEPTED MANUSCRIPT Stich, D., Serpelloni, E., Mancilla, F.d.L., Morales, J., 2006. Kinematics of the Iberia– Maghreb plate contact from seismic moment tensors and GPS observations. Tectonophysics 426, 295-317.
PT
Thurner, S., Palomeras, I., Levander, A.R., Carbonell, R., Lee, C.-T., 2014. Ongoing lithospheric removal in the western Mediterranean: Evidence from Ps receiver
SC RI
functions and thermobarometry of Neogene basalts (PICASSO project). Geochem. Geophys. Geosyst. 15, 1113-1127.
Timoulali, Y., Hahou, Y., Jabour, N., Merrouch, R., Kharrim, A., 2014. Main features of the deep structure by local earthquake tomography and active tectonics: case of
NU
Rif Mountain (morocco) and Betic Cordillera (Spain). J. Seismol. 18, 221-234. Torne, M., Fernandez, M., Comas, M.C., Soto, J.I., 2000. Lithospheric structure beneath
Geophys. Res. 105, 3209-3228.
MA
the Alboran basin: Results form 3D gravity modeling and tectonic relevand. J.
van Hinsbergen, D.J.J., Vissers, R.L.M., Spakman, W., 2014. Origin and consequences
ED
of western Mediterranean subduction, rollback, and slab segmentation. Tectonics 33, 393-419.
PT
Vergés, J., Fernández, M., 2012. Tethys–Atlantic interaction along the Iberia–Africa plate boundary: The Betic–Rif orogenic system. Tectonophysics 579, 144-172.
CE
Villasenor, A., Benz, H.M., Filippi, L., Scarpa, R., Patane, G., Vinciguerra, S., 1998. Three-dimensional P-wave velocity structure of Mt. Etna, Italy. Geophys. Res. Lett. 25, 1975-1978.
AC
Villaseñor, A., Gaite, B., Gallart, J., TOPO-IBERIA Seismic Working Group, 2011, Ambient noise tomography of the Iberian Peninsula using the IberArray seismic network, Geophysical Research Abstracts, Vol. 13, EGU2011-5853, EGU General Assembly 2011, Vienna.
23
ACCEPTED MANUSCRIPT
SC RI
PT
Figure captions
Figure 1. Map showing the structural units of the Iberia-Africa collision zone, and the seismic stations used in this study. B-R: Betic-Rif; Mz-d: Mesozoic deformed; Mz-u:
NU
Mesozoic undeformed. Thick lines indicate the deformation front of the Betic-Rif system and of the Middle Atlas (MA). BB: Beni Bousera peridotite Ce: Sarchal (Ceuta)
MA
peridotite; Gu: Guadalquivir basin; MM: Moroccan Meseta; Gh: Gharb Basin; RP: Ronda peridotite.
ED
Figure 2. Distribution of P-wave travel time residuals obtained using the IGN model as reference for two earthquakes (indicated by a purple star): a) a shallow depth event in
PT
northeastern Morocco and b) a subcrustal earthquake in the Gulf of Cadiz. Red circles
CE
indicate positive residuals (late/slow arrivals with respect to the 1-D IGN reference model) and blue circles indicate negative residuals (early/fast arrivals). Circle radius
AC
indicates the size of the residual according to the legend. Black lines indicate the boundaries between major structural units shown in Figure 1. IM: Iberian Massif; Gu: Guadalquivir Basin; EB: External Betics; IB: Internal Betics; IR; Internal Rif; ER: External Rif; MM: Moroccan Meseta; MA: Middle Atlas. Figure 3. Map of the model region and of the stations and earthquakes selected for traveltime tomography. Stations are indicated by open triangles and earthquakes by color-coded circles where color indicates focal depth according to the values in the palette. Al: Al-Hoceima; Ma: Málaga; Mo: Morón. The Trans-Alboran Shear Zone (TASZ) is indicated by the dashed line. The arrows show the left-lateral motion along the TASZ.
24
ACCEPTED MANUSCRIPT Figure 4. Indicators of the resolving power of the arrival time dataset used in this study for two layers of the model: 4-8 km (panels in left column) and 36-40 km (right
PT
column). a) Ray-path coverage (cumulative ray length through each cell in km) for the model layer from 4-8 km. Cells that are not sampled by any ray are shown in white. b)
SC RI
Same as a) but for the model layer from 36-40 km. c) Synthetic spike model for model layer 4-8 km. Velocity anomalies are shown in percentage with respect to the reference IGN model. Cells not sampled by the arrival time dataset are shown in gray. d) Same as
NU
c) but for the model layer from 36-40 km. e) Reconstruction of the synthetic model
MA
shown in c). f) Reconstruction of the synthetic model shown in d). Figure 5. Comparison of results before and after incorporation new arrival times from IberArray and other broadband stations. a) Horizontal slice of the model layer between
ED
8 and 12 km depth for the P-wave velocity model obtained using data from the IGN
PT
catalog alone for 2000-2009. White circles are earthquakes with focal depths inside this layer (after the simultaneous inversion for velocity structure and earthquake relocation).
CE
b) Same horizontal slice but for our preferred P-wave model obtained by combining
AC
IGN and IberArray arrival times. Note that the IGN-alone model only extends down to northernmost Morocco because the IGN bulletins include very few earthquakes farther south. P-wave velocity perturbations are shown in percentage with respect to the IGN 1D smoothed model. Cells not sampled by through-going rays are shown in grey. Figure 6. Horizontal slices of the obtained P-wave velocity model for different layers of the model. White circles are earthquakes with focal depths inside each layer (after the simultaneous inversion for velocity structure and earthquake relocation). Thin black lines indicate the boundaries between major structural units shown in Figures 1 and 2. P-wave velocity perturbations are shown in percentage with respect to the IGN 1-D smoothed model. Cells not sampled by through-going rays are shown in grey. a) Layer
25
ACCEPTED MANUSCRIPT 2 of the model, between 0 and 4 km depth. Thick line indicates the boundaries of the Granada-Guadix-Baza basins (Gr, Gu, B respectively). b) Layer 5 of the model,
PT
between 12-16 km depth. Thick lines indicate the location of the vertical cross sections
of the model, between 36 and 40 km depth.
SC RI
shown in Figure 7. c) Layer 9 of the model, between 28 and 32 km depth. d) Layer 11
Figure 7. Vertical cross sections of the P-wave velocity model. Earthquakes within 10 km of the cross-section are plotted as circles. The location of the cross sections is shown
NU
in Figure 6b. Vertical exaggeration is 2:1. a) NW-SE cross section through the Iberian
MA
Massif (IM), Ronda peridotite (Rp), Alboran basin (Alb) and Al-Hoceima (Al-H). b) WE cross section through the Gulf of Cadiz (GC) Strait of Gibraltar (SG), Alboran basin
AC
CE
PT
ED
(Alb), Alboran Ridge (AR) and western Algerian basin (Alg).
26
PT
ED
MA
NU
SC RI
PT
ACCEPTED MANUSCRIPT
AC
CE
Fig. 1
27
AC
CE
PT
ED
MA
NU
SC RI
PT
ACCEPTED MANUSCRIPT
Fig. 2
28
AC
Fig. 3
CE
PT
ED
MA
NU
SC RI
PT
ACCEPTED MANUSCRIPT
29
AC
CE
PT
ED
MA
NU
SC RI
PT
ACCEPTED MANUSCRIPT
Fig. 4
30
AC
CE
PT
ED
MA
NU
SC RI
PT
ACCEPTED MANUSCRIPT
Fig. 5
31
AC
Fig. 6
CE
PT
ED
MA
NU
SC RI
PT
ACCEPTED MANUSCRIPT
32
AC
CE
PT
ED
MA
NU
SC RI
PT
ACCEPTED MANUSCRIPT
Fig. 7
33
ACCEPTED MANUSCRIPT Table 1. IGN reference velocity model Limits (km)
Thickness (km)
1
0-11
11
2
11-24
13
3
24-31
halfspace
31-
AC
CE
PT
ED
MA
NU
SC RI
PT
Layer
34
Vp (km/s) 6.1 6.4
7
6.9
-
8.0
S0040-1951(15)00006-2 doi: 10.1016/j.tecto.2014.12.015 TECTO 126507
To appear in:
Tectonophysics
Received date: Revised date: Accepted date:
2 June 2014 14 November 2014 28 December 2014
Please cite this article as: El Moudnib, Lahcen, Villase˜ nor, Antonio, Harnafi, Mimoun, Gallart, Josep, Pazos, Antonio, Serrano, Inmaculada, C´ ordoba, Diego, Pulgar, Javier A., Ibarra, Pedro, Himmi, Mohammed M., Chourak, Mimoun, Crustal structure of the BeticRif system, western Mediterranean, from local earthquake tomography, Tectonophysics (2015), doi: 10.1016/j.tecto.2014.12.015
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT Crustal structure of the Betic-Rif system, western
SC RI
tomography
PT
Mediterranean, from local earthquake
Lahcen El Moudnib1,2 , Antonio Villaseñor3,* , Mimoun Harnafi1 , Josep
NU
Gallart3, Antonio Pazos4, Inmaculada Serrano5, Diego Córdoba6, Javier A.
MA
Pulgar7, Pedro Ibarra8, Mohammed M. Himmi2, and Mimoun Chourak9.
Institut Scientifique, Université Mohammed V-Agdal, Rabat, Morocco.
2
Faculté des Sciences, Université Mohammed V-Agdal, Rabat, Morocco.
3
Institute of Earth Sciences Jaume Almera, ICTJA-CSIC, Barcelona, Spain.
4
Real Instituto y Observatorio de la Armada, San Fernando, Cádiz, Spain.
5
Instituto Andaluz de Geofísica, Universidad de Granada, Granada, Spain.
6
IGEO and Universidad Complutense de Madrid, Madrid, Spain.
7
Departamento de Geodinámica Interna, Universidad de Oviedo, Oviedo, Spain.
8
Instituto Geológico y Minero de España, Madrid, Spain.
9
École Nationale des Sciences Appliquées, Université Mohammed Premier, Oujda,
AC
CE
PT
ED
1
Morocco.
*Corresponding author: [email protected]
Abstract We have determined the three-dimensional P-wave velocity structure of the crust and uppermost mantle beneath the Iberia-Africa collision zone using local earthquake
1
ACCEPTED MANUSCRIPT tomography. We have inverted arrival times of first-arriving P phases listed in the bulletins of the Instituto Geográfico Nacional and phases picked by us on permanent
PT
stations of other regional networks and on temporary broadband stations deployed in the frame of the TOPO-IBERIA project. In total we have used 38,624 P-wave arrival times
SC RI
from 2,362 local events recorded at 120 seismic stations. The most remarkable result is the imaging of a large high velocity body following approximately the western Alboran coastline, with a horizontal dimension of at least 200 km and extending in depth from
NU
the surface down to 24 km. This body, not imaged previously with this extent using
MA
seismic tomography, coincides with surface exposures of subcontinental mantle (peridotites) in Iberia and Africa and with a belt of positive gravity anomalies. We have also found a marked contrast in the seismic velocities of the middle and lower crust of
ED
the Alboran basin, coinciding with the location of the Trans-Alboran Shear Zone. We
PT
attribute this contrast to widespread magmatic intrusions in the eastern part of the basin, resulting in higher P-wave velocities than in the west. This contrast would also explain
CE
the different orientation of the Trans-Alboran Shear zone with respect to the surface
AC
features and faults in the Alboran basin. We also image thick crust beneath the Betics and Rif, accompanied by downgoing lithosphere of the Iberian foreland and Gulf of Cadiz beneath the Betic-Rif-Alboran system.
Keywords: tomography, Moho, Iberia, Morocco, Alboran
Highlights:
A high-velocity anomaly coincides with surface exposures of subcontinental mantle.
Different seismic structure between the western and eastern Alboran basin.
This velocity contrast might explain development of Trans-Alboran Shear Zone.
2
ACCEPTED MANUSCRIPT
The Iberian foreland basement deepens to the SE in the western Betics.
PT
1. Introduction The Iberia-Africa collision zone in the western Mediterranean is characterized by a
SC RI
complex interplate deformation zone. A tightly curved orogenic belt, the Betic-Rif system, is located between the undeformed Iberian and Morocco forelands (Figure 1). This section of the Eurasia-Africa plate boundary is characterized by a broad area of
NU
diffuse, moderate seismicity in contrast with the adjacent segments to the east (Algeria)
MA
and to the west (Gulf of Cadiz) where large earthquakes (M > 7) are frequent. The Betic-Rif system is traditionally divided into internal and external zones (e.g. Platt et al.,
ED
2013 and references therein). The external zones consist primarily of folded Mesozoic and Cenozoic sedimentary rocks from the former Iberian and African paleomargins,
PT
while the internal zones are composed of a stack of pre-Mesozoic and Mesozoic metamorphic rocks. Within the internal zones there are several exposures of lithospheric
CE
mantle with varying sizes: the Ronda peridotite in southern Spain is the largest outcrop
AC
of subcontinental mantle in the world with an extent up to 300 km2 (Obata, 1980); the Beni Bousera peridotite, with 70 km2 is the equivalent outcrop in the Rifean range; and finally some minor outcrops are also found in the Sarchal beach in Ceuta (SánchezGómez et al., 1995). The Alboran basin, located within the internal part of the Betic-Rif orogen, is an extensional basin formed during the Neogene (Comas et al., 1999). It is floored in the west by continental basement that largely corresponds to rocks belonging to the internal complexes of the Betics-Rif and by magmatic arc crust in the east (Booth-Rea et al., 2007). Although the geology and shallow structure of the region have been intensively studied (see recent reviews by Gutscher et al., 2012; Platt et al., 2013) the debate over the
3
ACCEPTED MANUSCRIPT kinematics and dynamics of its formation is still intense. Alternative evolution models based on kinematic reconstructions have been proposed for the region by Faccenna et
PT
al. (2004), Vergés and Fernández (2012) and van Hinsbergen et al. (2014) to mention only a few. Chertova et al. (2014) have tested these models using 4D thermomechanical
SC RI
numerical modeling and have been able to reproduce the slab morphology in the mantle observed using teleseismic tomography.
With respect to the deep structure of the crust and upper mantle, the first order features
NU
have been know for some time, for example the presence of a large high velocity
MA
anomaly in the upper mantle beneath the Alboran basin that probably corresponds to a subducted slab (Blanco and Spakman, 1993). Similarly, a significant number of studies have used local earthquake tomography with arrival times recorded by local networks to
ED
image the crustal structure of the region (Dañobeitia et al., 1998; Gurria and Mezcua,
PT
2000; Morales et al., 1999; Serrano et al., 2002; Serrano et al., 2003; Timoulali et al., 2014) . However, many of them are based on arrival times picked on analog recordings,
CE
and their resolution in the Alboran basin and northern Morocco is limited because of the
AC
sparse station coverage. In this study we make use of new high-quality data provided by the digitally-upgraded local seismic networks and also by a dense deployment of temporary broadband stations in southern Spain and northern Morocco. With this improved dataset we aim at obtaining more detailed images of the crust and uppermost mantle of the Iberia-Africa collision zone that can provide additional constraints on its past and present geodynamic evolution.
2. Data The data used in this study are first arrival times of P phases of local earthquakes, recorded at seismic stations. Because of the characteristics of the methodology, only
4
ACCEPTED MANUSCRIPT earthquakes and stations that are located inside the study region (model volume) are used. Our main source for arrival times is the bulletins of the Instituto Geográfico
PT
Nacional (IGN) from 2000-2009. The IGN operates the Spanish National Seismic Network (RSN, Red Sísmica Nacional), which by the year 2000 had already completed
SC RI
its upgrade to broadband sensors, digital recordings, and satellite communications. Therefore the phase arrival time data after 2000 are presumably of higher quality than earlier analog data. The IGN bulletins contain not only arrival times from stations of the
NU
Spanish RSN, but also include arrival times from neighboring networks such as the
MA
national Portuguese and Moroccan monitoring networks, operated by the Instituto Português do Mar e da Atmosfera (IPMA) and Centre National pour la Recherche Scientifique et Technique (CNRST) respectively. Additionally they include arrival time
ED
data from other research networks in the region, for example those operated by the
PT
Spanish Royal Navy Observatory (ROA). The station coverage provided by these networks is well suited for locating local earthquakes in the southern Iberian Peninsula,
CE
because of the greater station density in the most seismically active onshore areas.
AC
However this situation is less optimal for seismic tomography studies, for which a homogeneous station coverage would be preferable. The situation is still worse for northern Morocco because of the larger station spacing there and the lower quality of the arrival time data (for the time period considered here the national network was still analog). To ameliorate this situation the TOPO-IBERIA project, launched in 2007, included the deployment of the IberArray temporary network (Díaz et al., 2009; Díaz et al., 2010) with the objective of achieving a homogeneous station coverage in Spain with an interstation distance of 60 km or smaller. Simultaneously 18 temporary broadband stations with the same instrumentation as IberArray were also installed in northern
5
ACCEPTED MANUSCRIPT Morocco, in order to reduce the low station coverage provided by permanent networks. Subsequently other temporary networks (e.g. PICASSO, WILAS) have been installed in
PT
the Iberian Peninsula and Morocco both expanding and densifying the coverage of
SC RI
IberArray, but their data have not been incorporated in this study.
We have complemented the arrival time database of the IGN bulletins from 2000 to 2009 with phase picks made by us on digital recordings from the IberArray stations in
NU
Spain and Morocco for selected earthquakes in the time period 2007-2009. We have
MA
also incorporated picks made by us on recordings from other permanent networks that are not routinely included in the IGN bulletins, such as the Andalusian Seismic Network operated by the Instituto Andaluz de Geofísica (IAG). In total we have used data from
ED
120 stations in the region, 55 of them new broadband instruments previously
PT
unavailable (35 in the Iberian Peninsula and 20 in northern Africa). Figure 2 illustrates the improvement in coverage and quality provided by the new dataset. The use of a
CE
seismic array with a dense, homogeneous spacing not only improves the data coverage
AC
but also allows to identify errors, for example stations with residuals significantly different in size and/or sign than those of neighboring stations. In most cases these differences can be traced to bad pickings due to low signal-to-noise ratio or timing errors (stations that temporary lost the GPS time signal). The use of a seismic array also allows to evaluate qualitatively the strength and patterns present in the data. For example the distribution of residuals for the NE Morocco earthquake shown in Figure 2a clearly indicates a correlation with the different structural units. The Iberian and Moroccan forelands (Iberian Massif and Moroccan Meseta) and the Betics internal zones are clearly associated with negative (fast) residuals, while the external Betics and Rif exhibit large positive (slow) residuals. This general pattern is also observed for the
6
ACCEPTED MANUSCRIPT Gulf of Cadiz event shown in Figure 2b, although the negative residuals in the internal Betics and Portuguese Iberian Massif are now smaller in size or even positive. Stations
PT
located next to the Spanish Alboran shore are almost always consistently negative, indicating the presence of thinner than average crust and/or the presence of a large high-
SC RI
velocity body.
From this combined dataset we selected well recorded earthquakes for the tomographic
NU
inversion. First we choose a rectangular region of 880 720 km2 centered in the
MA
Alboran basin with a good distribution of local earthquakes and seismic stations (map in Figure 3 corresponds to the target region). Seismicity in this region is characterized by
ED
large variability in focal depth, including shallow, intermediate-depth and even a cluster of deep-focus earthquakes at approximately 650 km depth. Because of the decrease in
PT
seismicity with depth and in order to keep the model size manageable, we have considered only earthquakes with focal depth h < 100 km. We have then extracted
CE
earthquakes inside this volume with P phases recorded by 10 or more stations at
AC
distances less than 300 km, and with an azimuthal gap smaller than 200. Koulakov (2009) has shown that the use of events with azimuthal gaps greater than the commonly used value of 180 can contribute to improving the images obtained from local earthquake tomography. Here we use a value only slightly larger than 180 in order to try to improve the ray coverage in the marine parts of the model of interest, particularly in the Gulf of Cadiz. We have also required that selected earthquakes include 2 or more S readings because of their positive influence in improving the determination of focal depths (Gomberg et al., 1990). Finally, to avoid large outliers and misidentifications, we have eliminated earthquakes whose hypocenter solution reported by IGN had root mean square (RMS) residuals larger than 3.5 s. These selection criteria resulted in a dataset of
7
ACCEPTED MANUSCRIPT 38,624 P-wave arrival times from 2,362 local events recorded by 120 seismic stations, shown in Figure 3. Selected seismicity is rather well distributed over the entire model
PT
region, but three prominent seismic zones can be distinguished. The Trans-Alboran Shear Zone (TASZ) is a linear feature that extends approximately from Al-Hoceima in
SC RI
Morocco to southeastern Spain. A dense cluster of seismicity in the western Betics near the town of Morón is due to two moderate events (Mw ~ 4.5) in 2007-2008 and their aftershock sequences. While the seismicity in the region is mostly shallow (crustal), a
NU
remarkable zone of subcrustal, intermediate-depth seismicity occurs in the western
MA
Alboran basin, extending approximately N-S from Malaga to the west of Al-Hoceima. These deep earthquakes contribute to the dataset with ray-paths that improve the
ED
illumination of the uppermost mantle.
PT
3. Method
The tomographic inversion method used in this study is described in detail by Benz et
CE
al. (1996) and will only be summarized here. It has been successfully applied to
AC
different spatial scales, including volcanoes (Benz et al., 1996; Okubo et al., 1997; Villasenor et al., 1998) and to larger regions such as the one considered here (e.g. Powell et al., 2014). In this technique, arrival times of first-arriving P waves are inverted to simultaneously determine the three-dimensional P-wave velocity structure and earthquake relocations. This non-linear problem is linearized and solved iteratively. In each iteration, updated velocity model and earthquake locations are obtained and used as starting values for the next iteration. The iteration procedure is stopped after a preset number of iterations (in our case 10 iterations) or when the data misfit is not reduced significantly. Theoretical travel times are calculated using the finite-difference technique of (Podvin and Lecomte, 1991), which is well suited for computing travel
8
ACCEPTED MANUSCRIPT times in media with large lateral velocity perturbations. Velocity anomalies with respect to a reference model are solved using the iterative LSQR method of Paige and Saunders
PT
(1982) with smoothing constraints applied to control the degree of model roughness. We have tested a wide range of values of the smoothing parameter k (see Equation 7 in
SC RI
Benz et al. (1996) for its definition) in order to find those that produce good image fidelity without introducing large velocity alternations in adjacent cells and at the same time produce a significant reduction in data misfit. We have tested values of k ranging
NU
from 300 (smooth model) to 10 (very rough). Results show that decreasing the
MA
smoothing parameter down to a value of 75 continuously improves the data fit. Smaller values result in very rough models, and oscillations in the data fit with each iteration. Therefore we have adopted a value of k = 75, because it is the one that produces the best
ED
variance reduction without degrading the appearance of the model. A more detailed
PT
description of these tests is shown in the Supplementary Material, together with examples of under- and over-damped models.
CE
The reference/starting model considered here is the 1-D layered model used by the IGN
AC
for routine earthquake location in the Iberian Peninsula (Table 1). This model is representative of the velocity structure of the Iberian Massif (Paleozoic units in Figure 1) and has a crustal thickness of 31 km. We have interpolated the constant-velocity layers and discontinuities of the IGN 1-D model to linear gradients with smooth variations, in order to prevent that rays travel exclusively along the layer discontinuities and do not sample other parts of the model. However, if the arrival time data contains enough paths refracted through subsurface discontinuities (e.g. Pn phases) these discontinuities will naturally develop in the inverted tomographic model. The model volume, as described in the previous section, has horizontal dimensions of 880 720 km2 and, in order to include all potential ray-paths between earthquakes and
9
ACCEPTED MANUSCRIPT stations, it extends from 4 km above sea level to 120 km depth. The velocity model is parameterized as a regular grid of constant velocity cells of 20 20 4 km3 (in the x, y
PT
and z directions respectively), resulting in a total of 49,104 model parameters. A finer grid of 1 1 1 km3 was used for the finite-difference traveltime calculations. The size
SC RI
of the model and traveltime cells was determined empirically through a series of tests that addressed the effects of cell size and smoothing constraints (Benz et al., 1996).
NU
4. Results
MA
Inversion of the P-wave arrival time data after 10 iterations resulted in a reduction of 28% of the residual RMS from an starting value of 0.83 s. Because of the large number
ED
of unknowns, we cannot compute the formal model resolution and covariance, and therefore the resolution analysis of our model is done empirically using the ray-path
PT
coverage and synthetic reconstruction checks using the source-receiver geometry of our dataset. We use spike tests, that are similar to traditional checkerboard tests, except that
CE
the alternating positive and negative anomalies are not adjacent, but separated by a
AC
buffer with 0% velocity anomaly. These results show good recovery of the synthetic anomalies without significant horizontal smearing in the central part of the model, which is well sampled by through-going rays (Figure 4a,b). Resolution decreases below 40 km where most rays have vertical paths coming from intermediate-depth earthquakes in the western Alboran region. Hit count maps and results from spike tests for two layers of the model are shown in Figure 4 and with more detail in the Supplementary Material. As a result of the simultaneous inversion for velocity structure and earthquake location, new locations have been obtained. However, changes in location are not very large. The average difference in epicenter between the original IGN locations and the final
10
ACCEPTED MANUSCRIPT relocations is 4.5 km, and the average difference in depth is -1 km (the minus sign indicates that the final locations are shallower on average than the initial ones).
PT
In order to evaluate the influence of the added arrival times measured on IberArray stations in southern Spain and northern Morocco we have performed tomographic
SC RI
inversions with and without including this dataset. The result of this comparison is shown in Figure 5. The top panel shows a horizontal slice of the model (8-12 km depth) obtained using data listed in the IGN bulletins alone, while the bottom panel shows the
NU
same horizontal slice obtained using the combined dataset. An obvious improvement is
MA
the extension to the south of the region that is well covered by ray paths to include most of the Rif system in northern Morocco. There are, however some substantial differences even in regions that are reasonably well sampled by the IGN dataset. The most
ED
remarkable of these differences is in the shape and extension of the high velocity
PT
anomaly associated with the outcrop of the Ronda peridotites. While in the results using only the IGN bulletins this anomaly appears as a large blob, in the results from the
CE
combined dataset it has a curved, narrow, elongated shape, coincident with the Iberian
AC
coastline of the Alboran basin, which is in good agreement with the distribution of positive residual gravity anomalies in the region (Torne et al., 2000). Another relevant difference is the sharp contrast in velocity observed in the results of the combined dataset across the TASZ, with very low velocities in the west and faster velocities in the east, which is absent in the results based on the IGN bulletins alone. Figure 6 shows other horizontal layers of the model from the shallow crust to uppermost mantle. A common feature to all layers shown is the large variability in the values of the velocity perturbation with respect to the reference model, which often exceed 10%. This is particularly true for the negative (slow) velocity anomalies, which extend through large regions of the model.
11
ACCEPTED MANUSCRIPT Sedimentary basins are generally characterized by large low velocity anomalies. However, in our model the most important sedimentary basins (Guadalquivir, Gharb,
PT
and West Alboran basins) are poorly imaged because few stations and earthquakes are located on them, resulting in a lack of sampling rays. This is not the case for the
SC RI
intramontane basins in the Betic range, that are characterized by moderate seismicity and well monitored by seismic networks. The Granada-Guadix-Baza basins can be clearly identified as low velocities in the shallow layers of the model (Figure 6a).
NU
Similarly, the land portion of the Gulf of Cadiz accretionary complex is characterized
MA
by slow velocities, which are remarkably low in southernmost Spain, coinciding with the location of the Campo de Gibraltar flyschs. These low velocities not only occur at shallow depths, where they can be associated with the flysch units, but throughout the
ED
entire crust. Another conspicuous low velocity anomaly is located in the central Alboran
PT
basin and follows the trend of the TASZ. Although small sedimentary basins occur in this region, this anomaly is too large to be explained exclusively by sediments.
CE
The curved high velocity anomaly associated with the surface expression of the Ronda
AC
peridotite can be clearly followed from shallow depths (e.g. Figure 5b) to 12-16 km depth (Figure 6b) and down to at least 24 km (not shown here). This is also the case for the high velocity anomaly located east of the TASZ that is clearly imaged from 8 to 24 km depth. Figure 6c shows the layer of the model between 28 and 32 km depth, which is located just above the Moho of the reference model. A conspicuous feature in this layer is a high velocity zone that coincides with the extent of the Alboran basin. This is most likely caused by crustal thickness in this region being smaller than that of the reference model (31 km). Therefore the mantle at this depth would have higher velocities than the surrounding crustal blocks. The Ronda peridotite and eastern TASZ high velocity
12
ACCEPTED MANUSCRIPT anomalies described previously seem to be rooted in this high velocity anomaly,
PT
indicating that they extend throughout the crust down to the top of the mantle.
Figure 6d shows the layer between 36 and 40 km depth, within the mantle of the
SC RI
reference model. Large low velocity anomalies in this layer could indicate crustal velocities and therefore crustal thicknesses larger than the reference model. This is in fact observed all along the Betic-Rif system and would indicate the presence of a crustal
NU
root beneath those ranges.
MA
5. Discussion
ED
Here we will discuss the most relevant results obtained in this study and their
PT
implication for the geodynamic evolution of the region and present-day tectonics.
Anomalies associated with peridotite outcrops
CE
The most remarkable new finding in this study is the curved, narrow, high-velocity
AC
anomaly that follows the northwestern shore of the Alboran basin and coincides with the exposure of subcontinental mantle in the Ronda peridotite. Other studies based on local earthquake (Timoulali et al., 2014) and surface wave data (Chourak et al., 2005) have found high seismic velocities coinciding with the Ronda peridotites, but not with the geometry and geographical extent obtained here. This can be attributed to the lack of ray paths to stations in Morocco, and underscores the importance of the new dataset provided by the IberArray temporary stations. In our model this anomaly can be traced continuously for at least 200 km from approximately 3.5W in the east up to Ceuta in the southern shore of the Strait of Gibraltar. (e.g. Figure 6b), suggesting that the Ronda and Ceuta peridotite outcrops correspond to the same continuous body at depth.
13
ACCEPTED MANUSCRIPT Although the high velocity anomaly imaged in our model continues south of Ceuta it does not reach unequivocally the location of the Beni Bousera peridotites (Figure 6b).
PT
However other geophysical evidence such as gravity anomalies (Torne et al., 2000) and ambient noise tomography (Silveira et al., 2013; Villaseñor et al., 2011) suggest that
SC RI
this anomaly extends further south and therefore encompasses the three exposures of subcontinental mantle in the region. The inability of our model to image the continuation of the high velocity into the Beni Bousera peridotite is probably due to the
NU
lower station density in northern Morocco. As seen in Figure 1, few stations (none of
MA
them broadband) are located next to the Moroccan western Alboran shore nor in the Rif internal zones, preventing the detailed imaging of the structure in that region. Incorporation of other existing stations (such as those from the PICASSO experiment)
ED
would be necessary to improve the resolving power of the traveltime dataset in the
PT
internal Rif.
In the same way this high velocity anomaly has a large lateral extent, it also exhibits a
CE
large vertical continuity, extending from the surface to at least 24 km depth. This can be
AC
clearly observed in the vertical cross section of the model shown in Figure 7a. This figure also shows that the anomaly is not vertical but dips towards the center of the Alboran basin.
The Trans-Alboran Shear Zone (TASZ) Another new remarkable feature of the model is the strong change in P-wave velocity across the TASZ, with slow velocities in the west and fast in the east. This contrast, however is not observed at all depths, but most clearly between 8 and 24 km. The location of this boundary can be seen in the vertical cross section shown in Figure 7b. The orientation of the TASZ (N35E, Stich et al., 2006) is somewhat unexpected because
14
ACCEPTED MANUSCRIPT it is nearly perpendicular to the direction of convergence between Eurasia and Africa at this location, and is also different to the SW-NE orientation of the main sinistral faults
PT
in the Rif (Jebha and Nekor) and bathymetric features in the Alboran basin (e.g. the Alboran ridge). However, if the high velocity anomalies east of the TASZ correspond to
SC RI
a mechanically stronger block of the Alboran basin crust, it could act as an indenter causing the extrusion to the SW of a weaker western Alboran crust, thus explaining the left-lateral strike-slip mechanisms along the TASZ (see for example Figure 7 in Stich et
NU
al., 2003). Since the high velocity body is only observed below 8 km depth, it is
MA
compatible with structures in the upper crust having different orientations, for example inherited from the extension episode that generated the Alboran basin. The origin of this strong, high-velocity block, which coincides approximately with the eastern Alboran
ED
volcanic belt drawn by Duggen et al. (2003), could be magmatic intrusions in the
PT
preexisting thinned continental crust. In fact, magmatic arc crust is observed in the eastern Alboran basin, and arc-magmatism intrusions are more common in the eastern
CE
Alboran basin than in the west (Booth-Rea et al., 2007). The boundaries of crustal types
AC
determined by Booth-Rea et al. (2007) do not coincide exactly with the boundary of the block imaged here east of the TASZ but since those are based on shallow seismic studies, it is not incompatible with a different geometry at depth. The high velocity block seems to extend inland in Morocco (Figures 5b and 6b), in general agreement with previous tomographic studies (Serrano et al., 2003). In our model the boundary between high and low velocities seems to be located just east of Al-Hoceima, coinciding with an exposure of arc volcanic rocks. The inland extension of the western boundary of the high velocity block is also coincident with a region of a pronounced change in crustal thickness, with values of 25 km in the east to 40 km and deeper to the west (Mancilla et al., 2012).
15
ACCEPTED MANUSCRIPT
The onshore Gulf of Cadiz accretionary complex
PT
Very low velocities are observed at shallow depths in the westernmost external Betics. This anomaly coincides with the largest outcrop of flysch units in the region (e.g.
SC RI
Crespo-Blanc et al., 2012). Interestingly at all crustal depths of our model, this region is the one than exhibits the lowest seismic velocities. This observation has also been found independently by Silveira et al. (2013) using seismic ambient noise. In their study they
NU
obtain very low velocities (up to 40% slower than the average) for the propagation of
MA
Rayleigh waves travelling through the crust in the westernmost Betics and Rif. This indicates that P-wave velocities in the western Rif should be similarly slow to those in the western Betics if a comparable station and ray-path coverage were available.
ED
These anomalously low P- and surface-wave velocity values suggest that the flysch
PT
trough units and/or sediments of the Iberian-Morocco paleomargins have been underthrusted to great depths. This has been suggested for the Rif by Chalouan et al.
CE
(2001) and is also suggested by the cross sections in Figure 7 for the western Betics.
AC
Figure 7a clearly shows low-velocity rocks overridden by the high velocity anomaly associated to the Ronda peridotite. In map view there is also a good correlation between the location of this high velocity anomaly and the exposures of the flysch trough units (Figure 1) suggesting that the high velocity anomaly might have acted as a backstop for the underthrusting of the flyschs and other units.
The Iberian foreland basement Incorporation in the IGN bulletins of stations located in southern Portugal provide some constraints on the seismic structure of the southwestern Iberian Massif. Since the reference model used in this study (1-D model used by IGN for earthquake location) is
16
ACCEPTED MANUSCRIPT representative of the structure in the Iberian Massif, velocity anomalies in that region should be close to zero. This is indeed the case for the most of the crustal layers
PT
(Figures 6a, 5b and 6b). However, when we approach to the crust-mantle discontinuity some differences begin to emerge. In Figure 6c, corresponding to the layer of the model
SC RI
representing the reference Moho, we observe a substantial difference between the western and eastern Iberian Massif. Beneath Portugal the velocity anomalies continue to be close to zero or slightly negative. This behavior is already clearly present in the
NU
arrival time data (see the residual pattern in the Gulf of Cadiz earthquake shown in
MA
Figure 2b). Conversely, velocities in the eastern Iberian Massif are larger than average and extend farther south into the deformation front of the external Betics. It is not clear what is the reason for this east-west change in seismic velocities, because it does not
ED
coincide with any of the boundaries between the units that form the Iberian Massif.
PT
Crustal thickness throughout the Iberian Massif is remarkably constant, approximately 31 km (e.g. Díaz and Gallart, 2009), and this can be observed in the northern part of the
CE
cross section shown in Figure 7a. However, towards the SE there is a clear deepening of
AC
the contrast between crustal and mantle velocities. This deepening is also accompanied by a SE dipping zone of intermediate-depth earthquakes, that can be followed up to the surface. Similarly, although not as well constrained as for the Iberian Massif, the crustmantle discontinuity beneath the Gulf of Cadiz and the earthquake focal depths seem to dip towards the east (Figure 7b). This suggest that the continental lithosphere of Iberia (and also the lithosphere beneath the Gulf of Cadiz) is undergoing continental subduction (or alternatively delamination or roll-back) beneath the Alboran Sea. Furthermore, in this part of the western Betics, the Iberian lithosphere would still be attached to the slab that has been imaged in the mantle beneath the Alboran Sea (e.g. Garcia-Castellanos and Villaseñor, 2011; Spakman and Wortel, 2004). A similar
17
ACCEPTED MANUSCRIPT observation has been obtained using receiver functions (Mancilla et al., 2013). This process of continental subduction/delamination could also be present in the Rif, but
SC RI
Crustal thickness beneath the Betics and Rif
PT
cannot be constrained with our present dataset.
Regions in which the crust is thicker than the reference Moho would be characterized in Figure 6c by significant low velocity anomalies, which correspond exactly to the entire
NU
Betic-Rif system, including areas such as the Strait of Gibraltar that are not
MA
characterized by high topography. A very similar pattern in crustal thickness is also observed using receiver functions (Thurner et al., 2014; Mancilla et al., submitted), and ambient noise tomography (Villaseñor et al., 2011). Between 36-40 km depth (Figure
ED
6d) strong low velocities are still observed beneath the central Betics, indicating the
CE
6. Conclusions
PT
existence of a deeper root beneath the areas of highest topography.
AC
We have obtained a new P-wave velocity model of the Iberia-Africa collision zone, that includes the Betic-Rif system and the Alboran basin. A major finding of this study is a large high velocity body 200-km long and extending from 0 to approximately 24 km depth that is associated with the surface exposures of subcontinental mantle in the Ronda and Sarchal (Ceuta) peridotites, and possibly extending to the Beni Bousera peridotites in the Rif. Because of its size and location, this body, not previously imaged with this extent by similar studies, might hold important clues with respect to the evolution of the Alboran-Betic-Rif region. Likewise, because of its potential effect in teleseismic travel times and surface wave dispersion, it should be taken into account for seismic studies of the mantle beneath this region that do not consider crustal structure
18
ACCEPTED MANUSCRIPT and use instead crustal corrections. We have also imaged the different character of the eastern and western Alboran middle and lower crust, and its potential implications for
PT
the development of the TASZ. Most of these results have been possible because of the incorporation of arrival times
SC RI
measured on the temporary IberArray broadband stations in southern Spain and Morocco. Due to the lower station density in northern Morocco (although greatly improved with respect to previous studies) some features such as the continuation of the
NU
high velocity body into the Beni Bousera peridotite, the estimate of crustal thickness
MA
beneath the Rif, and the behavior of the African foreland basement towards the Alboran basin have not been completely resolved. Incorporation of data from other existing temporary deployments in the region (both marine and on land) might help to resolve
PT
ED
these issues.
CE
Acknowledgements
(IGN, Spain) for providing the bulletin and
AC
waveform data that are the base of this study. We also thank other networks operating in the region that have provided waveform data for this study: Instituto Andaluz de Geofísica (IAG, Spain), Real Instituto y Observatorio de la Armada (ROA, Spain), Instituto Português do Mar e da Atmosfera (IPMA, Portugal). Data from other individual stations have been obtained from the ORFEUS (The Netherlands), GEOFON (Germany) and IRIS (USA) data centers. This is a contribution of the Team ConsoliderIngenio 2010 TOPO-IBERIA (CSD2006-00041). Additional funding was provided by the SIBERIA (CGL2006-01171), RIFSIS (CGL2009-09727) and ALERTES-RIM (CGL2013-45724-C3-3-R) projects.
19
ACCEPTED MANUSCRIPT References Benz, H.M., Chouet, B.A., Dawson, P., Lahr, J.C., Page, R.A., Hole, J.A., 1996. Threedimensional P and S wave velocity structure of Redoubt Volcano, Alaska J.
PT
Geophys. Res. 101, 8111-8128.
Blanco, M.J., Spakman, W., 1993. The P-wave velocity structure of the mantle below
Tectonophysics 221, 13-34. -
-
SC RI
the Iberian Peninsula: evidence for subducted lithosphere below southern Spain.
, J.M., Grevemeyer, I., 2007. Crustal
types and Tertiary tectonic evolution of the Alborán sea, western Mediterranean.
NU
Geochem. Geophys. Geosyst. 8, Q10005.
Chalouan, A., Michard, A., Feinberg, H., Montigny, R., Saddiqi, O., 2001. The Rif
MA
mountain building (Morocco): a new tectonic scenario. Bulletin de la Société Géologique de France 172, 603-616.
Chertova, M.V., Spakman, W., Geenen, T., van den Berg, A.P., van Hinsbergen, D.J.J.,
ED
2014. Underpinning tectonic reconstructions of the western Mediterranean region with dynamic slab evolution from 3-D numerical modeling. Journal of Geophysical
PT
Research: Solid Earth 119, 5876-5902. Chourak, M., Corchete, V., Badal, J., Gómez, F., Serón, J., 2005. Shallow Seismic
CE
Velocity Structure of the Betic Cordillera (Southern Spain) from Modelling of Rayleigh Wave Dispersion. Surv. Geophys. 26, 481-504.
AC
Comas, M.C., Platt, J.P., Soto, J.I., Watts, A.B., 1999. The origin and tectonic history of the Alboran Basin: insights from LEG 161 results, in: Zahn, R., Comas, M.C., Klaus, A. (Eds.). Ocean Drilling Program, College Station, TX, USA, pp. 555-580. Crespo-Blanc, A., Balanyá, J.C., Expósito, I., Luján, M., Suades, E., 2012. Crescentlike large-scale structures in the external zones of the western Gibraltar Arc (BeticRif orogenic wedge). Journal of the Geological Society 169, 667-679. Dañobeitia, J.J., Sallarès, V., Gallart, J., 1998. Local earthquakes seismic tomography in the Betic Cordillera (southern Spain). Earth Planet. Sci. Lett. 160, 225-239. Díaz, J., Gallart, J., 2009. Crustal structure beneath the Iberian Peninsula and surrounding waters: A new compilation of deep seismic sounding results. Phys. Earth Planet. Inter. 173, 181-190. Díaz, J., Villaseñor, A., Gallart, J., Morales, J., Pazos, A., Córdoba, D., Pulgar, J.A., García-Lobón, J.L., Harnafi, M., Group, T.-I.S.W., 2009. The IBERARRAY
20
ACCEPTED MANUSCRIPT broadband seismic network: a new tool to investigate the deep structure beneath Iberia. Orfeus Newsl. 8, 1-6. Díaz, J., Villaseñor, A., Morales, J., Pazos, A., Córdoba, D., Pulgar, J.A., Garcia-Lobon,
PT
J.L., Harnafi, M., Carbonell, R., Gallart, J., 2010. Background Noise Characteristics at the IberArray Broadband Seismic Network. Bull. Seism. Soc. Am. 100, 618-628.
SC RI
Duggen, S., Hoernle, K., van den Bogaard, P., Rüpke, L.H., Phipps Morgan, J., 2003. Deep roots of the Messinian salinity crisis. Nature 422, 602-606. Faccenna, C., Piromallo, C., Crespo-Blanc, A., Jolivet, L., Rosseti, F., 2004. Lateral slab deformation and the origin of the western Mediterranean arcs. Tectonics 23,
NU
TC1012.
Garcia-Castellanos, D., Villaseñor, A., 2011. Messinian salinity crisis regulated by
MA
competing tectonics and erosion at the Gibraltar arc. Nature 480, 359-363. Gomberg, J.S., Shedlock, K.M., Roecker, S.W., 1990. The effect of S-wave arrival times on the accuracy of hypocenter estimation. Bull. Seism. Soc. Am. 80, 1605-
ED
1628.
Gurria, E., Mezcua, J., 2000. Seismic tomography of the crust and lithospheric mantle
PT
in the Betic Cordillera and Alboran Sea. Tectonophysics 329, 99-119. Gutscher, M.A., Dominguez, S., Westbrook, G.K., Le Roy, P., Rosas, F., Duarte, J.C.,
CE
Terrinha, P., Miranda, J.M., Graindorge, D., Gailler, A., Sallares, V., Bartolome, R., 2012. The Gibraltar subduction: A decade of new geophysical data. Tectonophysics 574-575, 72-91.
AC
Koulakov, I., 2009. Out-of-Network Events Can Be of Great Importance for Improving Results of Local Earthquake Tomography. Bull. Seism. Soc. Am. 99, 2556-2563. Mancilla, F., Stich, D., Berrocoso, M., Martin, R., Morales, J., Fernández-Ros, A., Paez, R., Pérez-Peña, A., 2013. Delamination in the Betic Range: Deep structure, seismicity, and GPS motion. Geology 41, 307-310. Mancilla, F.d.L., Stich, D., Morales, J., Julià, J., Díaz, J., Pazos, A., Córdoba, D., Pulgar, J.A., Ibarra, P., Harnafi, M., Gonzalez-Lodeiro, F., 2012. Crustal thickness variations in northern Morocco. J. Geophys. Res. 117, B02312. Morales, J., Serrano, I., Jabaloy, A., Galindo-Zaldívar, J., Zhao, D., Torcal, F., Vidal, F., González-Lodeiro, F., 1999. Active continental subduction beneath the Betic Cordillera and the Albor
. Geology 27, 735-738.
21
ACCEPTED MANUSCRIPT Obata, M., 1980. The Ronda Peridotite: Garnet-, Spinel-, and Plagioclase-Lherzolite Facies and the P-T Trajectories of a High-Temperature Mantle Intrusion. J. Petrol. 21, 533-572.
PT
Okubo, P.G., Benz, H.M., Chouet, B.A., 1997. Imaging the crustal magma sources beneath Mauna Loa and Kilauea volcanoes, Hawaii. Geology 25, 867-870.
SC RI
Paige, C.C., Saunders, M.A., 1982. LSQR: An algorithm for sparse linear equations and sparse least squares. ACM Trans. Math. Software 8, 43-71.
Platt, J.P., Behr, W.M., Johanesen, K., Williams, J.R., 2013. The Betic-Rif Arc and Its Orogenic Hinterland: A Review. Annu. Rev. Earth Planet. Sci. 41, 313-357.
NU
Podvin, P., Lecomte, I., 1991. Finite difference computation of traveltimes in very contrasted velocity models: a massively parallel approach and its associated tools.
MA
Geophys. J. Int. 105, 271-284.
Powell, C.A., Withers, M.M., Cox, R.T., Vlahovic, G., Arroucau, P., 2014. Crustal velocity structure associated with the eastern Tennessee seismic zone: Vp and Vs
ED
images based upon local earthquake tomography. Journal of Geophysical Research: Solid Earth 119, 464-489.
PT
Sánchez-Gómez, M., García Dueñas, V., Muñoz, M., Balanyá, J.C., 1995. Relación estructural entre los cuerpos peridotíticos situados al Norte y al Sur del Estrecho de
CE
Gibraltar. Geogaceta 17, 135-137. Serrano, I., Zhao, D., Morales, J., 2002. 3-D crustal structure of the extensional Granada Basin in the convergent boundary between the Eurasian and African plates.
AC
Tectonophysics 344, 61-79. Serrano, I., Zhao, D., Morales, J., Torcal, F., 2003. Seismic tomography from local crustal earthquakes beneath eastern Rif Mountains of Morocco. Tectonophysics 367, 187-201. Silveira, G., Dias, N.A., Villaseñor, A., 2013. Seismic imaging of the western Iberian crust using ambient noise: Boundaries and internal structure of the Iberian Massif. Tectonophysics 589, 186-194. Spakman, W., Wortel, M.J.R., 2004. A tomographic view on Western Mediterranean Geodynamics, in: Cavazza, W., Roure, F., Spakman, W., Stampfli, G.M., Ziegler, P.A. (Eds.), The TRANSMED Atlas, The Mediterranean Region from Crust to Mantle, pp. 31-52. Stich, D., Ammon, C.J., Morales, J., 2003. Moment tensor solutions for small and moderate earthquakes in the Ibero-Maghreb region. J. Geophys. Res. 108, 2148. 22
ACCEPTED MANUSCRIPT Stich, D., Serpelloni, E., Mancilla, F.d.L., Morales, J., 2006. Kinematics of the Iberia– Maghreb plate contact from seismic moment tensors and GPS observations. Tectonophysics 426, 295-317.
PT
Thurner, S., Palomeras, I., Levander, A.R., Carbonell, R., Lee, C.-T., 2014. Ongoing lithospheric removal in the western Mediterranean: Evidence from Ps receiver
SC RI
functions and thermobarometry of Neogene basalts (PICASSO project). Geochem. Geophys. Geosyst. 15, 1113-1127.
Timoulali, Y., Hahou, Y., Jabour, N., Merrouch, R., Kharrim, A., 2014. Main features of the deep structure by local earthquake tomography and active tectonics: case of
NU
Rif Mountain (morocco) and Betic Cordillera (Spain). J. Seismol. 18, 221-234. Torne, M., Fernandez, M., Comas, M.C., Soto, J.I., 2000. Lithospheric structure beneath
Geophys. Res. 105, 3209-3228.
MA
the Alboran basin: Results form 3D gravity modeling and tectonic relevand. J.
van Hinsbergen, D.J.J., Vissers, R.L.M., Spakman, W., 2014. Origin and consequences
ED
of western Mediterranean subduction, rollback, and slab segmentation. Tectonics 33, 393-419.
PT
Vergés, J., Fernández, M., 2012. Tethys–Atlantic interaction along the Iberia–Africa plate boundary: The Betic–Rif orogenic system. Tectonophysics 579, 144-172.
CE
Villasenor, A., Benz, H.M., Filippi, L., Scarpa, R., Patane, G., Vinciguerra, S., 1998. Three-dimensional P-wave velocity structure of Mt. Etna, Italy. Geophys. Res. Lett. 25, 1975-1978.
AC
Villaseñor, A., Gaite, B., Gallart, J., TOPO-IBERIA Seismic Working Group, 2011, Ambient noise tomography of the Iberian Peninsula using the IberArray seismic network, Geophysical Research Abstracts, Vol. 13, EGU2011-5853, EGU General Assembly 2011, Vienna.
23
ACCEPTED MANUSCRIPT
SC RI
PT
Figure captions
Figure 1. Map showing the structural units of the Iberia-Africa collision zone, and the seismic stations used in this study. B-R: Betic-Rif; Mz-d: Mesozoic deformed; Mz-u:
NU
Mesozoic undeformed. Thick lines indicate the deformation front of the Betic-Rif system and of the Middle Atlas (MA). BB: Beni Bousera peridotite Ce: Sarchal (Ceuta)
MA
peridotite; Gu: Guadalquivir basin; MM: Moroccan Meseta; Gh: Gharb Basin; RP: Ronda peridotite.
ED
Figure 2. Distribution of P-wave travel time residuals obtained using the IGN model as reference for two earthquakes (indicated by a purple star): a) a shallow depth event in
PT
northeastern Morocco and b) a subcrustal earthquake in the Gulf of Cadiz. Red circles
CE
indicate positive residuals (late/slow arrivals with respect to the 1-D IGN reference model) and blue circles indicate negative residuals (early/fast arrivals). Circle radius
AC
indicates the size of the residual according to the legend. Black lines indicate the boundaries between major structural units shown in Figure 1. IM: Iberian Massif; Gu: Guadalquivir Basin; EB: External Betics; IB: Internal Betics; IR; Internal Rif; ER: External Rif; MM: Moroccan Meseta; MA: Middle Atlas. Figure 3. Map of the model region and of the stations and earthquakes selected for traveltime tomography. Stations are indicated by open triangles and earthquakes by color-coded circles where color indicates focal depth according to the values in the palette. Al: Al-Hoceima; Ma: Málaga; Mo: Morón. The Trans-Alboran Shear Zone (TASZ) is indicated by the dashed line. The arrows show the left-lateral motion along the TASZ.
24
ACCEPTED MANUSCRIPT Figure 4. Indicators of the resolving power of the arrival time dataset used in this study for two layers of the model: 4-8 km (panels in left column) and 36-40 km (right
PT
column). a) Ray-path coverage (cumulative ray length through each cell in km) for the model layer from 4-8 km. Cells that are not sampled by any ray are shown in white. b)
SC RI
Same as a) but for the model layer from 36-40 km. c) Synthetic spike model for model layer 4-8 km. Velocity anomalies are shown in percentage with respect to the reference IGN model. Cells not sampled by the arrival time dataset are shown in gray. d) Same as
NU
c) but for the model layer from 36-40 km. e) Reconstruction of the synthetic model
MA
shown in c). f) Reconstruction of the synthetic model shown in d). Figure 5. Comparison of results before and after incorporation new arrival times from IberArray and other broadband stations. a) Horizontal slice of the model layer between
ED
8 and 12 km depth for the P-wave velocity model obtained using data from the IGN
PT
catalog alone for 2000-2009. White circles are earthquakes with focal depths inside this layer (after the simultaneous inversion for velocity structure and earthquake relocation).
CE
b) Same horizontal slice but for our preferred P-wave model obtained by combining
AC
IGN and IberArray arrival times. Note that the IGN-alone model only extends down to northernmost Morocco because the IGN bulletins include very few earthquakes farther south. P-wave velocity perturbations are shown in percentage with respect to the IGN 1D smoothed model. Cells not sampled by through-going rays are shown in grey. Figure 6. Horizontal slices of the obtained P-wave velocity model for different layers of the model. White circles are earthquakes with focal depths inside each layer (after the simultaneous inversion for velocity structure and earthquake relocation). Thin black lines indicate the boundaries between major structural units shown in Figures 1 and 2. P-wave velocity perturbations are shown in percentage with respect to the IGN 1-D smoothed model. Cells not sampled by through-going rays are shown in grey. a) Layer
25
ACCEPTED MANUSCRIPT 2 of the model, between 0 and 4 km depth. Thick line indicates the boundaries of the Granada-Guadix-Baza basins (Gr, Gu, B respectively). b) Layer 5 of the model,
PT
between 12-16 km depth. Thick lines indicate the location of the vertical cross sections
of the model, between 36 and 40 km depth.
SC RI
shown in Figure 7. c) Layer 9 of the model, between 28 and 32 km depth. d) Layer 11
Figure 7. Vertical cross sections of the P-wave velocity model. Earthquakes within 10 km of the cross-section are plotted as circles. The location of the cross sections is shown
NU
in Figure 6b. Vertical exaggeration is 2:1. a) NW-SE cross section through the Iberian
MA
Massif (IM), Ronda peridotite (Rp), Alboran basin (Alb) and Al-Hoceima (Al-H). b) WE cross section through the Gulf of Cadiz (GC) Strait of Gibraltar (SG), Alboran basin
AC
CE
PT
ED
(Alb), Alboran Ridge (AR) and western Algerian basin (Alg).
26
PT
ED
MA
NU
SC RI
PT
ACCEPTED MANUSCRIPT
AC
CE
Fig. 1
27
AC
CE
PT
ED
MA
NU
SC RI
PT
ACCEPTED MANUSCRIPT
Fig. 2
28
AC
Fig. 3
CE
PT
ED
MA
NU
SC RI
PT
ACCEPTED MANUSCRIPT
29
AC
CE
PT
ED
MA
NU
SC RI
PT
ACCEPTED MANUSCRIPT
Fig. 4
30
AC
CE
PT
ED
MA
NU
SC RI
PT
ACCEPTED MANUSCRIPT
Fig. 5
31
AC
Fig. 6
CE
PT
ED
MA
NU
SC RI
PT
ACCEPTED MANUSCRIPT
32
AC
CE
PT
ED
MA
NU
SC RI
PT
ACCEPTED MANUSCRIPT
Fig. 7
33
ACCEPTED MANUSCRIPT Table 1. IGN reference velocity model Limits (km)
Thickness (km)
1
0-11
11
2
11-24
13
3
24-31
halfspace
31-
AC
CE
PT
ED
MA
NU
SC RI
PT
Layer
34
Vp (km/s) 6.1 6.4
7
6.9
-
8.0