RESEARCH ARTICLE
10.1002/2016GC006657
Lithospheric structure of Iberia and Morocco using
finite-frequency Rayleigh wave tomography from earthquakes
and seismic ambient noise
I. Palomeras
1,2
, A. Villase
~
nor
3
, S. Thurner
1
, A. Levander
1
, J. Gallart
3
, and M. Harnafi
4
1
Rice University, Houston, Texas, USA,
2
Department of Geology, University of Salamanca, Salamanca, Spain,
3
Institute of
Earth Sciences Jaume Almera, ICTJA-CSIC, Barcelona, Spain,
4
Scientific Institute of Rabat, Universit
e Mohammed V-Agdal,
Rabat, Morocco
Abstract We present a new 3-D shear velocity model of the western Mediterranean from the Pyrenees,
Spain, to the Atlas Mountains, Morocco, and the estimated crustal and lithospheric thickness. The velocity
model shows different crustal and lithospheric velocities for the Variscan provinces, those which have been
affected by Alpine deformation, and those which are actively deforming. The Iberian Massif has detectable
differences in crustal thickness that can be related to the evolution of the Variscan orogen in Iberia. Areas
affected by Alpine deformation have generally lower velocities in the upper and lower crust than the Iberian
Massif. Beneath the Gibraltar Strait and surrounding areas, the crustal thickness is greater than 50 km, below
which a high-velocity anomaly (>4.5 km/s) is mapped to depths greater than 200 km. We identify this as a
subducted remnant of the NeoTethys plate referred to as the Alboran and western Mediterranean slab.
Beneath the adjacent Betic and Rif Mountains, the Alboran slab is still attached to the base of the crust,
depressing it, and ultimately delaminating the lower crust and mantle lithosphere as the slab sinks. Under
the adjacent continents, the Alboran slab is surrounded by low upper mantle shear wave velocities
(Vs < 4.3) that we interpret as asthenosphere that has replaced the continental margin lithosphere which
was viscously removed by Alboran plate subduction. The southernmost part of the model features an
anomalously thin lithosphere beneath the Atlas Mountains that could be related to lateral flow induced by
the Alboran slab.
1. Introduction and Tectonic Setting
The westernmost Mediterranean comprises the Iberian Peninsula and Morocco, separated by the Alboran
Sea and the Algerian Basin (Figure 1). This area, the far western end of the Alpine-Himalayan orogenic belt,
is currently affected by the Africa-Eurasia convergence with deformation extending from the Pyrenees in
the north of Iberia to the Atlas Mountains in Morocco. The current convergence rate is 3 mm/yr with Iberia
moving to the southwest relative to Africa [e.g., Koulali et al., 2011].
A series of tectonic events have affected this area since the Paleozoic. From late Ordovician to the Permian,
different continental blocks were accreted to form the supercontinent Pangaea. The Variscan orogeny
resulted from the collision of Baltica-Laurentia and Gondwana in Devonian to Carboniferous times [Matte,
1986]. During continental convergence, the Armorica and Avalonia microplates docked against Baltica-
Laurentia in the active margin with Gondwana [Matte, 2001]. In Iberia, this orogeny is represented by the
Iberian Massif which outcrops in almost the entire western half of the peninsula (Figure 1a). In Africa, Varis-
can age Gondwanan rocks are found in the Moroccan Meseta, the Anti-Atlas, and the Sahara craton. During
the Mesozoic, Pangaea broke up as the Neo-Tethys Ocean opened from east to west, separating the future
Eurasian plate from Africa. Breakup included formation of rift grabens and basins in both northern Africa
and Iberia. Iberia became a microplate surrounded by shallow waters and sedimentary basins. In the Middle
Jurassic to early Cretaceous, the Atlantic opened to the west of Iberia, with the Iberian plate rotating
counter-clockwise and opening the Bay of Biscay. Initiation of African-Eurasian convergence in the Creta-
ceous accreted Iberia to the Eurasian plate forming the Pyrenees. Continued African-Eurasian convergence
is thought to have transmitted the plate margin stress fields that reached the interior of the Iberian plate,
forming a thrust system on the previous Mesozoic rift basins, creating, e.g., the Iberian Chain. The
Key Points:
3-D Shear velocity model of Iberia
and north Morocco
Moho and LAB map of Iberia and
north Morocco
Crustal and Upper mantle imaging
Correspondence to:
I. Palomeras,
imma@usal.es
Citation:
Palomeras, I., A. Villase
~
nor, S. Thurner,
A. Levander, J. Gallart, and M. Harnafi
(2017), Lithospheric structure of Iberia
and Morocco using finite-frequency
Rayleigh wave tomogra phy from
earthquakes and seismic ambient
noise, Geochem. Geophys. Geosyst., 18,
1824–1840, doi:10.1002/
2016GC006657.
Received 20 SEP 2016
Accepted 13 MAR 2017
Accepted article online 28 MAR 2017
Published online 4 MAY 2017
V
C
2017. American Geophysical Union.
All Rights Reserved.
PALOMERAS ET AL. RAYLEIGH WAVE TOMOGRAPHY IN WM 1824
Geochemistry, Geophysics, Geosystems
PUBLICATIONS
topography of the Central System and that of the Cantabrian Mountains together with their foreland basins,
the Duero and Tajo tertiary basins were also formed during this time. In North Africa, a Mesozoic Rift System
was inverted and uplifted since the Middle Eocene, forming the Atlas Mountains [Jacobshagen et al., 1988;
Laville et al., 2004]. Cenozoic African-European convergence resulted in subduction of the Tethys oceanic
plate beneath parts of southern Europe. The subduction zone against northeastern Spain began to retreat
rapidly starting at 30–25 Ma, separating into multiple segments as it traversed eastward toward Italy and
southward to Africa, forming the modern western Mediterranean basins [Royden, 1993; Rosenbaum et al.,
2002; Chertova et al., 2014; Van Hinsbergen et al., 2014]. One of these segments moved southwestward,
opening the Algerian and Alboran basins. The retreating trench reached the passive margins of south Iberia
and north Morocco in the Middle-Late Miocene, initiating construction of the Betics-Rif fold and thrust belts
and their foreland basins, the Guadalquivir in Iberia and the Rharb in Morocco [Rosenbaum et al., 2002; Van
Hinsbergen et al., 2014].
Neogene volcanic fields are found in northeastern Iberia (Catalan Volcanic Province, CVP), the Valencia
Trough, the eastern Betics and Rif, the Alboran Sea, the isolated Calatrava Volcanic Field (CVF) north of the
Betics, and are widespread in the Middle and High Atlas. All these volcanics, except the CVF and the ones in
the Atlas, show a similar evolution from silica rich magmas in the Early to Middle Miocene, to silica poor
magmas from the Late Miocene to present [e.g., Mart
ıetal., 1992; Duggen et al., 2005]. The change in com-
position indicates a shift from calc-alkaline subduction-related volcanism to asthenosphere sourced melts.
The geochemistry of the CVF and the Atlas volcanic provinces indicates their alkaline magmas have a sub-
continental lithospheric origin [e.g., Cebri
a and L
opez-Ruiz, 1995; Bosch et al., 2014].
Figure 1. (a) Tectonic map of the Western Mediterranean. Orange triangles represent the Cenozoic volcanic fields. The Iberian Massif terranes have been named: SPZ: South Portuguese
Zone; OMZ: Ossa-Morena Zone; CIZ: Central Iberian Zone; WALZ: West Asturian-Leonese Zone; CZ: Cantabrian Zone. (b) Topographic map with the broadband stations used in this study.
Color codes are used to signify different arrays (inset).
Geochemistry, Geophysics, Geosystems 10.1002/2016GC006657
PALOMERAS ET AL. RAYLEIGH WAVE TOMOGRAPHY IN WM 1825
Following the terminology of Gibbons and Moreno [2002] for Spain, Iberia, and Morocco can be divided into
(1) Variscan zones consisting of the largely stable Iberian Massif and Moroccan Meseta, where Paleozoic
rocks outcrop and (2) the Cenozoic Alpine zones of eastern and southern Iberia, and those from Morocco
(the Rif and Atlas Mountains) (Figure 1a). Currently the western Mediterranean is under compressive defor-
mation due to the convergence between Eurasia and Africa. Active deformation is distributed over the
Betic-Rif chain, Gulf of Cadiz, and Alboran Sea as a diffuse band [Buforn and Udias, 2010] with high seismic
activity in the crust and at intermediate depths, and a few deep events at 650 km depth beneath Granada
[e.g., Buforn et al., 1991]. Local tectonics in this region are complicated with extension occurring in the east-
ern Betics and Rif, compression to the west, and strike slip on the Trans Alboran Shear Zone (TASZ) [Stich
et al., 2006]. The active region is surrounded by the less deformed Variscan zone in Iberia and by the Sahara
Craton to the south of the Atlas Mountains.
Body wave tomography studies have imaged a high-velocity anomaly beneath the westernmost Mediterra-
nean that has been interpreted as the remnant of trench rollback of the Alboran slab [ Blanco and Spakman,
1993; Calvert et al., 2000; Gutscher et al., 2002; Garcia-Castellanos and Villase
~
nor, 2011; Bezada et al., 2013].
Recent active and passive seismic studies of the lithosphere have identified unusually thick crust under the
western Betics and Rif [Gil et al., 2014; Palomeras et al., 2014; Thurner et al., 2014] and imaged crustal delami-
nation beneath the central Betics and Rif [Palomeras et al., 2014; Thurner et al., 2014] in response to the
Alboran slab subduction. Lithosphere thickness variations suggest that the subducting Alboran slab viscous-
ly thins the adjacent continental margin lithosphere mantle [Levander et al., 2014].
Our goal in this paper is to present an improved lithosphere and asthenosphere shear velocity model of the
western Mediterranean from the Pyrenees to the Sahara craton, using Rayleigh wave phase velocity disper-
sion. We used the extensive data set available from recent experiments and permanent networks. The avail-
ability of this dense array data provides a unique opportunity to obtain high-resolution images of the
western Mediterranean lithosphere.
2. Data and Methodology
We used data from 368 broadband seismic stations of permanent and temporary arrays deployed in the
area (Figure 1b). Permanent networks included those of Morocco, Spain, and Portugal. Temporary arrays
included the recent IberArray (2007–2013, Spain), PICASSO (2009–2012, USA), M
€
unster University (2010–
2012, Germany), and Bristol University (2010–2012, UK) experiments. Average interstation spacing is approx-
imately 60 km (Figure 1b) with almost uniform coverage across Iberia and northern Morocco, but less regu-
lar coverage in the Middle and High Atlas Mountains. To obtain a velocity model from the surface to
200 km depth, we did Rayleigh wave tomography using both ambient noise and teleseismic earthquakes.
2.1. Ambient Noise Tomography
We have obtained phase velocity maps of fundamental mode Rayleigh waves between 4 and 40 s from
cross correlation of seismic ambient noise. To obtain the maps presented here we followed the approach
used by Silveira et al. [2013]. The main difference with respect to these authors is that we have used a longer
time period to compute the stacked cross correlations (2008–2011), and a much larger data set that includ-
ed 350 of the broadband stations described previously.
The first step for ambient noise tomography is the determination of empirical Green’s functions by cross
correlating and stacking continuous recordings of all the potential station pairs. We follow the procedure of
Bensen et al. [2007]: first, we extracted 4 h time windows from the continuous recording, eliminating those
that contained 2 or more data gaps, filled the remaining gaps with the average of the signal before and
after the gap, converted all the records to velocity, and decimated them to five samples per second. In order
to avoid the contamination of the noise recordings by earthquakes and nonstationary sources near the sta-
tions, we have applied a ‘‘temporal normalization’’ by dividing the noise signal by the running average (with
a width of 256 s) of the absolute value of the waveform amplitude filtered between 20 and 100 s (the domi-
nant period band of the earthquake surface wave energy). Then we have applied spectral normalization or
‘‘whitening’’ to correct for the fact that the ambient noise amplitude spectrum is not flat in the potential
period band of interest (e.g., 1–50 s). This is achieved by dividing the true amplitude spectrum of each noise
window by the smoothed amplitude spectrum. Finally, with the objective of improving the signal-to-noise
Geochemistry, Geophysics, Geosystems 10.1002/2016GC006657
PALOMERAS ET AL. RAYLEIGH WAVE TOMOGRAPHY IN WM 1826
ratio, all the processed recordings are cross correlated and stacked for the entire common time period avail-
able to each station pair.
The empirical Green’s functions constructed using this methodology contain causal and acausal Rayleigh
waves (in positive and negative times, respectively). To measure phase velocities we calculated the symmet-
ric cross correlation by averaging the causal and acausal segments (reversing the time for the latter). Phase
velocity dispersion is then measured using an automated implementation of the frequency-time analysis
(FTAN) methodology [Levshin et al., 1992]. Following Bensen et al. [2007], from the complete set of measure-
ments we have selected for tomography those phase velocity measurements that have values of signal-to-
noise ratio (SNR) > 10.
Finally, we applied a finite-frequency tomographic inversion [Barmin et al., 2001] to the selected phase
velocity measurements between 4 and 40 s periods to obtain fundamental-mode Rayleigh-wave phase
velocity maps on a 0.58 3 0.58 grid across the study area. The tomography method of Barmin et al. [2001]
approximately accounts for the spatially extended frequency-dependent sensitivity of the waves by using
Gaussian sensitivity kernels.
2.2. Teleseismic Surface-Wave Tomography
We analyzed 168 teleseismic events occurring between April 2009 and December 2011 with magnitudes
greater than 6.0 and epicentral distances between 308 and 1208 that were recorded by a large fraction of
the 368 available broadband stations. Not all the stations recorded the entire time period as the IberArray
stations were relocated as part of the operational plan [D
ıaz et al., 2009]. The path density is high at all peri-
ods (Figure 2b), and the azimuthal distribution of events is good, with the largest number of events coming
from the northeast and the lowest from the SSE (Figure 2a).
We analyzed the vertical component Rayleigh wave signals. The seismograms were corrected for instrument
response, and filtered with 18 Butterworth band-pass filters in the period band of 20–167 s (20, 22, 25, 27,
30, 34, 40, 45, 50, 59, 67, 77, 87, 100, 111, 125, 143, and 167 s). The fundamental mode Rayleigh wave is iso-
lated from other modes and body waves by windowing filtered records with a variable length tapered win-
dow. The amplitude was corrected for frequency-dependent anelastic attenuation and geometrical
spreading [Mitchell, 1995] and normalized to the root mean square amplitude for each event.
The 2-D phase velocity for periods from 20 to 167 s was calculated using the modified two-plane wave tech-
nique as described by Forsyth and Li [2005] with 2-D finite-frequency kernels for both amplitude and phase
[Yang and Forsyth, 2006]. To preserve the plane wave assumption, the study area was divided into two
regions. The northern region comprises the Iberian Peninsula while the southern one comprises the Betics,
the Alboran Sea, and the Atlas mountains (Figure 2b). The phase velocities of the southern region are
described in Palomeras et al. [2014]. We followed the same procedure and parameters to obtain the phase
velocity for the northern region. Although the difference between the calculated phase velocities in the
overlapping portions (from 2108 to 0.58 in longitude and from 358 to 428 in latitude) is not greater than
0.10 km/s in most of the periods and points, decreasing to less than 0.05 km/s to larger periods, the phase
velocities were merged using a weighted averaging function. At the northern edge of the southern box, the
weight is zero, and increases linearly to one moving south. The same weighting is applied to the northern
region increasing from zero to one as we move north. Both areas have the same weight value (0.5) at the
central point of the overlap.
3. Phase Velocity Maps
The phase velocity maps (Figure 3) from both the ANT and the teleseismic Rayleigh waves show higher val-
ues in Variscan Iberia than in Alpine Iberia at short periods (4–35 s and 25–40 s, respectively). These results
image the consequences that the age of tectonic events and the orogenic evolution of the Iberian Peninsu-
la have in its crustal structure. We observe unusually low-phase velocities beneath Gibraltar from 4 to 50 s
suggesting a thick crust. Low-phase velocities are also found in the eastern end of the Rif and the Betics,
and beneath the Middle and High Atlas in the period range 40 s up to 87 s. These periods typically sample
the lowermost crust and uppermost mantle.
Geochemistry, Geophysics, Geosystems 10.1002/2016GC006657
PALOMERAS ET AL. RAYLEIGH WAVE TOMOGRAPHY IN WM 1827
4. Inversion
The ANT and teleseismic phase velocities data sets were calculated on a regular 0.58 3 0.58 grid and overlap
in periods from 20 to 40 s. To obtain the final dispersion curves from 4 to 167 s, we averaged the phase
velocity for the overlapping periods (25–35 s) at each grid node. Following Yao et al. [2008], we averaged
the phase velocities, and found all overlapping phase velocities differed by less than 0.15 km/s. Actually, the
80% of the phase velocities differ less than 0.10 km/s. At 20 s, we used the values obtained from ANT and at
Figure 2. (a) Azimuthal distribution of earthquakes (red dots) used in this study. (b) Ray coverage shown as a gray scale for 50 s period in the two regions into which the study area has
been divided (see text). Open red triangles are the stations used in this study. White dashed lines outline the areas with good ray coverage, those with at least 1400 hits.
Geochemistry, Geophysics, Geosystems 10.1002/2016GC006657
PALOMERAS ET AL. RAYLEIGH WAVE TOMOGRAPHY IN WM 1828
