The P-wave velocity structure of the mantle below the Iberian Peninsula: evidence for subducted lithosphere below southern Spain

Tectonophysics, 221 (1993) 13-34

13

Elsevier Science Publishers B.V., Amsterdam

The P-wave velocity structure of the mantle below the Iberian Peninsula: evidence for subducted lithosphere below southern Spain Maria JosC Blanc0

a and Wim Spakman

b

a Departamento

de Geofuica, Uniuersidad Complutense de Madrid, Spain b Department of Geophysics, University of Utrecht, The Netherlands

(Received March 13, 1991; revised version accepted September 16, 1991)

ABSTRACT Blanco, M.J. and Spakman, W., 1993. The P-wave velocity structure of the mantle below the Iberian Peninsula: evidence for subducted lithosphere below southern Spain. In: J. Badal, J. Gallart and H. Paulssen (Editors), Seismic Studies of the Iberian Peninsula. Tectonophysics, 221: 13-34. The P-wave velocity structure of the mantle below the Iberian Peninsula is investigated with the method of delay-time tomography. More than 200,000 delay-times are inverted for simultaneous estimates of P-wave velocity anomalies, event mislocation and origin time errors, and station corrections. Below the Betic-Alboran area we find a positive anomaly between 200 and 700 km depth which in shape and amplitude resembles the image of a subducted lithosphere slab. Laterally this anomaly exhibits a clear SW-NE trend. In view of the low level of spatial resolution encountered we focus in our analysis on assessing the reliability of the slab-like image. Arguments in favor of the existence of a slab structure in the upper mantle are developed from inversions of subsets of the data. The resolving power for imaging a slab below the Betic-Alboran region is investigated using sensitivity analyses with a velocity model for a subducted slab. We conclude that a slab-structure exists in the mantle below the Betic and Alboran region and that it can be imaged with sufficient reliability. We interpret the anomaly as the image of subducted lithosphere. The slab is presumably detached from the surface. We propose that subduction took place during at least part of the Oligocene followed by slab detachment in the early Miocene.

Introduction Most of the seismological studies of the Iberian Peninsula (Fig. 1) have resulted in information about the crust and the lithosphere-asthenosphere system (e.g., Payo, 1970, 1972; Banda et al., 1980, 1981, 1983; Panza et al., 1980; Marillier and Mueller, 1985; Snieder, 1988; Badal et al., 1990; Paulssen and Visser, 1993-this issue). Some insight in the structure of the deeper mantle is obtained from seismic tomography studies (Spakman, 1986, 1990, 1991; Granet and Trampert, 1989; Nolet, 1990). The tomographic results ob-

Correspondence to: Dr. M.J. Blanco, Instituto Geografico Nacional, Avd. Anaga, 35; planta 11, 38001 Santa Cruz de Tenerife, Tenerife, Spain. Fax: 34 22 243017.

0040-1951/93/$06.00

tained below the Iberian Peninsula have only a limited resolution for structural detail, i.e., detail on the scale of 200-300 km or less seems poorly resolved, which hampers the interpretation of the mantle images in terms of mantle processes. Still, it is in the mantle where we may find important leads for studying the large scale evolution of the Iberian Peninsula. For instance, evidence for active deep mantle processes is given by the peculiar deep seismicity at about 640 km (e.g., Buforn et al., 1991). For lack of alternatives, the cause of these events is generally ascribed to the dynamics of a detached piece of cold lithosphere (e.g., Chung and Kanamori, 1976). It is called detached because there is no mantle seismicity between a depth of about 150 km and the depth of the deep events. A different indication that a subducted slab exists below Iberia was found in a tomo-

0 1993 - Elsevier Science Publishers B.V. All rights reserved

14

M.J. BLANC0

Data selection and delay-time tomography We apply the method of delay-time tomography developed by Spakman and Nolet (1988) and Spakman (1988, 1991) to study the P-wave velocity structure of the mantle beneath the Iberian Peninsula. P delay-time data were selected from the ISC bulletin tapes (years 1964-19861, but only from events that are recorded regionally or glob-

Peninsula

and surrounding

to a depth

of 1420 km. The epicenters

The panel

to the right shows the location

from events

and stations

located

regions.

displayed

outside

Below the area displayed,

in the Left panel

of regional

stations

the region

displayed.

W. SPAKMAN

ally by at least 10 seismological stations. Regional events contributed data from stations up to 90” of epicentral distance, i.e., from both regional and distant stations. For distant events, i.e., with epicenters outside the boundaries of the region of interest (Fig. 11, we only selected the delay-times observed in the regional stations. Over 200,000 delay-times were selected within these limits. Notice that we combine data from regional and teleseismic events in our analysis. The celf model employed to parameterize the slowness field of the mantle consists of 1” X 1” cells arranged in 20 layers to a depth of 1420 km. Each layer consists of 18 X 18 cells. The thickness of the cells increases gradually with depth from 33 km for the crust and lithosphere to 100 km in the lower mantle to a depth of 1420 km (see Table 1). The delay-times are inverted for simultaneous estimates of (1) cell slowness anomalies, (2) event mislocation and origin time errors, and (3) station delay-time corrections. For teleseismic events, only one time correction is computed instead of all four event parameters. We remark that the

graphic study by Spakman (1990). The slab-image obtained is, however, of questionable accuracy due to poor sampling of the Iberian mantle by the ray path data used, and also because the eastern part of the Iberian Peninsula was located near the edge of the investigated area. These results and the puzzling dynamics in the deeper mantle motivated us to perform a more detailed tomographic investigation of the P velocity structure of the mantle below the Iberian Peninsula, using data from both regional and teleseismic events. Specifically, we are interested in locating the mantle structure which is related to the deep seismicity.

Fig. 1. The lberian

AND

belong

the mantle velocity

structure

to the local subset of the earthquakes

is investigated

from which data are used. Note that also delay-time Solid lines within units.

the Iberian

Peninsula

down

used in this study.

delineate

data are used major

tectonic

P-WAVE VELOCITY STRUCTURE

OF THE MANTLE: EVIDENCE FOR SUBDUCTED

TABLE I Cell layer division and average velocity of model PM2 km

km

km/s

0 33 70 120 170 220 275 330 390 460 530 600 670 740 820 920 1020 1120 1220 1320

33 70 120 170 220 27.5 330 390 460 530 600 670 740 820 920 1020 1120 1220 1320 1420

6.228 7,853 8.028 8.117 8.109 8.176 8.414 8.732 9.319 9.710 9.941 10.151 10.899 11.016 11.183 11.360 Il.535 11.686 11.841 11.988

event parameters pertain to clusters of events and are, hence, average estimates for the mislocation and origin time error of a number of events near to each other (see Spakman and Nolet, 1988, for this approach). In our case, we discriminate between regional and teleseismic event clusters. Regional event clusters are defined as all events in cells of a size of 0.5” X 0.5” X 35 km, where “regional” is used here to denote all events within 30” of distance to the center of the Iberian Peninsula. To define teleseismic event clusters, we used larger cells of a size of 2.5” x 2.5” x 100 km. Of course, these cells stand in no relation to the cell model used for tomography; they are merely used to identify event clusters. In the inversion the only distinction between regional and teleseismic events lies in the number of unknowns to parameterize the event mislocation (see also Spakman and Nolet, 1988, and Spakman, 1991, for how we combine event parameters and velocity heterogeneity in a tomographic inversion). Because we invert for event and station parameters, there is no need to correct the data prior to inversion with, e.g., the average event delays and/or the average station delays. Hence,

15

LITHOSPHERE BELOW S SPAIN

we work with the actual ISC delay-times and not with so-called relative residuals. Only one correction of the data is necessarily applied prior to inversion and is related to the particular choice of the reference velocity model. The reference model is used in tomography to linearize the model equations, and for the computation of the ray path geometry, and the relocation and station correction coefficients (e.g., Spakman, 1991). The accuracy of the reference model, in the sense of approximating the laterally averaged Earth structure, is of utmost importance for the reliability of the tomographic results (Van der Hilst and Spakman, 1989; Zielhuis et al., 1989). The approach developed by the latter authors was used by Spakman et al. (1991) to derive-from 2.1 million ISC delays-a reference model, called PM2 (Fig. 21, valid for P delay-time tomography of the Eu-

6

P-velocity (km/s) 7

8

9

10

11

12

Fig. 2. Radial P-velocity models: Jeffrey+Bullen (JB) velocity and PM2. The latter model is determined by Spakman et al. (1991) and serves as a better reference model for delay-time tomography of the European-Mediterranean mantle including the mantle below Iberian Peninsula. See aiso Table 1.

16

M.J. BLANC0

ropean-Mediterranean mantle including the mantle below the Iberian Peninsula. In our study we will use model PM2 Since the ISC data have been computed using the Jeffreys-Bullen travel time tables, we have to correct the data prior to inversion for the difference in travel-time computed from model PM2 and the Jeffreys-Bullen tables. To compute this correction we still use the ISC event location, because any mislocation of the (ISO events in the PM2 model can be accounted for by the event parameters included in the inversion (Van der Hilst, 1990). A density-plot of the corrected data as a function of epicentral distance and delay-time value is shown in Figure 3. The data are more or less evenly distributed about the line of zero delay which corresponds to the predicted reference travel time in model PM2. The inversion of the tomographic equations is performed with the LSQR algorithm of Paige

delay ti~~s:96000

AND

W. SPAKMAN

and Saunders (19821, introduced in seismic tomography by Nolet (1985). We follow the inversion procedure as described by Spakman and Nolet (1988). The inversions are slightly damped only to make sure that the smallest eigenvalues of the problem do not have large effect on the solution, and lateral smoothing of cell slowness amplitudes is applied during inversion. In addition to delay-time inversions, we also performed sensitivity analyses with a number of synthetic velocity models. The synthetic tests are carried out to obtain estimates of accuracy and spatial resolution. Inversion results The first inversion we perform is of all the delay-time data with values between -3 and +3 s. We take these delay-time limits to exclude data

1

128

Fig. 3. Density-plot of delay-time value vs. epicentral distance. Delay-time counts are made in “boxes” of 0.1 s by 0.5”. The contouring is logarithmic. The delay-time data are computed relative to model PM2. The line of zero delay corresponds to the reference travel time in model PM2. Only data between - 3 and + 3 s (lines) are used in this study.

P-WAVE

VELOCITY

STRUCTURE

OF THE

MANTLE:

EVIDENCE

FOR

which possibly pertain to triplication branches or to exclude data that are possibly contaminated by large data errors. A data fit (reduction in Rit4S value of the difference between the actual delaytimes and those predicted from the 3-D inversion result) of 20% is obtained in 25 iterations with the LSQR algorithm. We will address this tomographic solution as result IPl; it is presented in Figure 4. The P-wave velocity anomalies are contoured in percentage deviations from the average cell layer velocities in reference model PM2 (Table 1). For judging the significance of the mapped anomalies we have to study the accuracy of the mapping. To this purpose we have computed a cell-spike sensitivity result. The cell-spike model consists of equidistantly spaced single-cell anomalies with a constant amplitude of 5% velocity contrast with the ambient reference mantle velocity (Table 1). Synthetic data are computed from the spike model along the same ray paths as used for the data inversions. We have added normally distributed noise with a standard deviation of 0.6 s to the synthetic data prior to inversion. This choice of the noise level resulted in a data fit of about 20%; similar to the data fit obtained for result IPl. We note that in the synthetic data inversion the event relocation and station correction parameters are also included in the model equations. The cell-spike inversion result is displayed in Figure 5. The regular grid of shaded circles indicates the actual locations and spatial size of the single-cell anomalies. The sensitivity test result is line-contoured in 0.5% anomaly value steps (i.e., in 10% steps of the peak value). The synthetic velocity model is only well retrieved in parts of the cell layers which are centered at depths of 51,95 and 145 km. In other parts of the model, specifically in the deeper layers, at most 20% of the spike value seems recoverable and the resolution for detail is poor or even absent. From this we conclude that the resolving power of the ray path set for detailed imaging of anomalies, i.e., with an accuracy on the order of the cell size, is generally poor and we need not consider interpreting small scale anomaly variations on the order of the cell size in IPl, except in the layers centered at 51, 95 and 145 km depth.

SUBDUCI-ED

LITHOSPHERE

BELOW

S SPAIN

17

The poor ability to map anomaly detail does not exclude the possibility that anomaly patterns of any distance scale can be detected. The poor resolution indicates that the mapping is heavily blurred in the sense that cell-scale anomalies are not detected or are smeared over many cells and may even be displaced relative to their actual locations, and that larger scale anomalies may be detected but are likely of unreliable shape and also smeared over many cells in, as yet, unknown directions, Therefore, result IPl may still bear important information on the structure of the mantle beneath the Iberian Peninsula which we have to retrieve by further analysis. To study the resolution for larger detail we perform a second kind of spike test in which the synthetic model consists of large blocks (2 x 2 x 2 adjacent cells) again with a synthetic amplitude of 5%. Normally distributed noise with a standard deviation of 0.8 s is added to the synthetic data in order to mimic the convergence of the real data inversion IPl. The result of this test is shown in Figure 6. We infer that the resolution for detail on the order of 200-300 km is considerably better. At many places single blocks can be identified in the inversion result. Still, the resolution is far from ideal considering the underestimated amplitudes and the anomaly smearing between many blocks. Nevertheless, this test demonstrates that the resolving power of the ray set employed may be sufficient to detect, in places resolve, mantle anomalies in the order of 200-300 km or larger. Hence, we collected sufficient evidence to motivate further study of the large scale patterns in IPl. An interesting anomaly pattern in IPl, which we will further investigate, is the positive anomaly located partly under the Alboran Sea and the Betic region at depths between 200 and 700 km and trending SW-NE. Our interpretation of this anomaly is that of a blurred (i.e., poorly resolved) image of a subducted slab. Subducted slabs of relatively aged lithosphere are mapped in tomography with positive anomaly amplitudes which primarily can be attributed to the temperature difference between slab and ambient mantle. Of course, before establishing such an interpretation we have to investigate the existence and reliability of the “slab anomaly”. A first indication that

M.J. BLANC0

1Data inversion: IPI Fig. 4. Tomographic

result IPl: velocity

1

anomalies

+2*0%

-2.0%

relative

AND W. SPAKMAN

to model PM2. Numbers

indicate

the depth at the center

of a cell-layer.

P-WAVE

VELOCITY

STRUCTURE

OF THE MANTLE:

EVIDENCE

FOR SUBDUCTED

LI-lXOSPHERE

BtLOW

S SPAIN

19

lSpike inversion: IPlJ Fig. 5. Spike sensltwlty contouring

result.

Grey

spots

is in 0.5% anomaly

denote

the actual

value increments.

location Dashed

of the cell-spikes. contour

lines indicate

The spike

anomaly

value

negative

anomaly

values.

is -i-5%. The

M.J.

20

Block inversion Fig. 6. Block sensitivity

result.

AND

W. SPAKMAN

1

Shaded

value is + 5%. The contouring

BLANC0

areas denote

is in 0.5% anomaly

the actual

location

value increments.

and spatial Dashed

size of the block anomalies.

contour

lines indicate

negative

The block anomaly anomaly

values.

P-WAVE

VELOCITY

STRUCTURE

OF THE

MANTLE:

EVIDENCE

FOR

the resolution may be sufficient to map a subducted slab below the Betic-Alboran region stems from the block inversion result. In Figure 7 we display two cross-sections through the block model which demonstrate that single blocks are reasonably resolved in depth and that not much anomaly smearing occurs across the 670 km boundary between upper and lower mantle. On itself this may be taken as sufficient indication that the “slab anomaly” is a resolvable, at least detectable, feature. However, since the sensitivi~ tests only give an indication about the resolving power of the ray set employed and not of the actual delay-time data (and actual ray paths) it is interesting to try to derive the existence of this anomaly (i.e., existence below the Betic-Alboran region and in the upper mantle) from the actual data. Furthermore, it is interesting to design a sensitivity test with a synthetic slab model to

SUBDUCTED

LITHOSPHERE

BELOW

S SPAIN

21

demonstrate explicitly to what extent a subducted slab can be resolved with our reference-model ray set. A similar test has been used by Spakman et al. (1989) to study the resolution for subducting slabs in the NW Pacific. To these purposes we will try to answer the following two questions: (1) is the “slab anomaly” a mapping of heterogeneity located in the upper mantle below the Betic and Alboran Sea’?, and (2) can a subducted slab located in the upper mantle below the Betic and Alboran Sea be imaged with sufficient reliability? Locating the “slab anomaly” To answer the first question, we performed additional inversions with two different data sets. The first data set is restricted to teleseismic data (48-90”) and is thus confined to combinations of either regional stations and distant events or re-

200

Fig. 7. Two cross-sections through the block sensitivity result of Figure 6. Shaded areas denote the actual location and spatial size of the block anomalies. The block anomaly value is + 5%. The contouring is in 0.25% anomaly value increments. Dashed contour lines indicate negative anomaly values.

M.J. BLANC0

[Data inversion: Fig. 8. Tomographic

IP2

1

result IP2 determined epicentral

W. SPAKMAN

+2.0%

-2.0%

from inverting distances

AND

the delay times from regional

between

48 and 90” (“teleseismic

station data”).

and distant

event combinations

with

P-WAVE.

VELOCITY

STRUCTURE

OF THE

MANTLE:

EVIDENCE

(Data inversion: IP31 Fig. 9. Tomographic

resutt

IP3 determined

FOR

SUBIXJCTED

LITHOSPHERE

BELOW

S SPAIN

23

-2.0% from

distances

inverting

the delay

times

less than 30” (“short-distance

from

station-event

data”).

combinations

with

epicentral

24

gional events and distant stations. The ray paths pertaining to this data set sample the upper mantle of the source (station) region, the lower mantle between source (station) region and the cell model, and the rays enter (leave) the cell model in its lower mantle extension. The inversion result is called IP2 and is displayed in Figure 8. The second data set consists of the data from station-event combinations with epicentral distance smaller than 30”. The ray paths belonging to these delay-time data sample the upper 750 km of the mantle in and outside the cell model. The inversion result will be addressed as result IP3 and is shown in Figure 9. Except for the upper mantle volume within the cell model, the two ray paths sets sample different mantle regions. The data gap of 30-48” ensures this (the gap of 18” corresponds to the lateral dimensions of the cell model). We will use the characteristics of the ray sets to demonstrate that the “slab anomaly” has an origin in the upper mantle. We remark that the resolving power of the data sets used to compute IPl, IP2, and IP3 is different. From sensitivity analyses with cellspikes similar to that shown for IPl in Figure 5 we have inferred that the IP2 result is slightly better resolved in the lower mantle than IPl and poorer in the upper mantle, and that IP3 is slightly better resolved in the upper part of the upper mantle and poorer in the lower part of the upper mantle. The resolution of IPl can be seen as a combination of these results; the resolution is slightly impaired in the uppermost mantle and the lower mantle, but slightly improved in the central part of the cell model. However, differences in resolution are small. In both results IP2 and IP3 we find a positive anomaly below the Betic and Alboran Sea. The “slab anomalies” in IPl and IP2 compare best. In IP3 the anomaly has a somewhat different shape but is still detected. Let us consider the “slab anomaly” in some detail. In Figure 10 we have plotted 2 cross-sections through the core of the anomaly for inversion result IPl, IP2, and IP3. One section (left) is taken perpendicular to the apparent strike direction of the anomaly. The other section runs parallel to the strike. The “slab anomaly” in results IPl and IP2 is rather

M.J. BLANC0

AND

W. SPAKMAN

similar. In IP3 the anomaly has an irregular shape and consists of two parts which, as we will discuss below, is due to lack of resolving power of the IP3 data set in the lower part of the upper mantle. In general, we find that mapped anomaly patterns exhibit somewhat elongated shapes. This apparent elongation is a result of poor resolution which causes smearing of anomalies in the main directions of ray sampling. This can be exemplified by studying the cell hit-count (Fig. 11). The cell hit-count is the number of rays sampling a particular cell. The three-dimensional variation in the ray sampling influences the inversion result strongly (see Spakman and Nolet, 1988; Spakman, 1988). The main directions of ray-illumination can be inferred from the main trends in the hit count and these are in-line with the direction of apparent elongation of the anomalies. One could demonstrate these trends explicitly by computing the ray density tensor GCissling, 1988) which contains all directional information. For our purposes the hit-count variation as displayed in Figure 11 suffices. Specifically, we can use this information to explain why the “slab anomaly” in IP3 deviates from that in IP2. Since IP3 is constructed from short distance data only, the ray paths sample the upper mantle at relative shallow angles (see hit-count pattern). The “slab anomaly” is smeared along these shallow angle directions, quite opposite to the smearing in IP2. Furthermore, the lower part of the “slab anomaly” is poorly illuminated by rays (low hit-count) which affects the resolution even more, i.e., in negative sense. These inferences are examples of the general statement we made in the previous section with regard to the poor resolution. The rays from which IP2 is constructed enter the cell model in the lower mantle. The rays from which IP3 is constructed only sample the upper mantle. The cross-sectional volume of both raysets is the upper mantle inside the cell model. As the “slab anomaly” is present in both results IP2 and IP3, we have found an important indication that the structural cause of the “slab anomaly” is at least located in the upper mantle of the cell model, a lower mantle extension is not excluded. But, one could ask: could the slab anomaly be a

P-WAVh

VELOCI’I’Y

STRUCTURE

0

OF THE

500

MANTLE:

EVIDENCE

FOR SUBDUCTED

1500

1000

2000

0

LITHOSPHERk2

500

BEl,OW

2s

S SPAIN

1000

1500

2000

400 g c g

600 800 foO0 1200

200 400 5 r g

600 800 1000 1200

200 400 g

600 800 1000 1200

Fig. 10. Two cross-sections each

column

apparent

of panels

strike

through indicates

the tomographic the exact

of the slab anomaly.

indicate

the locations

results

location

The other

of earthquakes.

IPl, 1P2, and IP3. The horizontal

of the section.

section

runs along

The cross-section strike

straight

line in the map at the top of

to the left is taken

of the anomaly.

Dots

Notice the deep event at 640 km near the center

perpendicular

in the maps

to the

and cross-sections

of the cross sections.

M.J. BLANC0

500

0

1500

1000

AND

W. SPAKMAN

2000

200 400 g z

600

B

800 1000 1200

200 400 E

600

f 0 D

800 1000 1200

Fig.

11. The cell hit-count

hit-count

of

10,000

ray-sampling

rays.

for the same Note

that

sections

the

of the cells. For reference

cell

as in Figure

hit-count

we plotted

10. The contouring

displayed the contours

for

the

sections

of the positive

scale

is logarithmic,

is computed anomalies

from

of Figure

i.e. 4 corresponds the

to a

three-dimensional

10 with thick lines.

P-WAVE

VELOCITY

STRUCTURE

OF THE MANTLE:

EVIDENCE

FOR SUBDUCTED

mapping of velocity structures outside the cell model or could it be a projection of station delay-time averages into cell anomalies, all due to poor resolution. Our answer is that this is not very likely. Recall that we also estimate in the inversions the station delay-time corrections and event mislocation parameters for all stations and events. For those stations and events outside the cell model these model parameters also absorb, in a way, the effects of mantle heterogeneity between the cell model and station location, and between the cell model and event, respectively (Spakman and Nolet, 19881. This considerably reduces the likelyhood that the “‘slab anomaly” is an erroneous mapping of velocity structure actually located outside the cell model or a projection of station averages on to the slowness cells. Moreover, any artifacts related to velocity structures located outside the cell model would rather be confined to those cells which are located in the vicinity of the edge of the cell model instead of creating a large amplitude anomaly in the upper central part of the cell model. The average station delays for stations in the Betic and along the eastern margin of Spain are invariably positive (Fig. 12) indicating that negative velocity anomalies would be expected beneath the stations (if the station residuals can be related to velocity heterogeneity in the first place), and that an erroneous mapping of the station delay into celI-anomalies will not result into an anomaly with a positive amplitude but into one with a negative amplitude. The arguments put forward above strongly support the existence of a high-velocity body in the upper mantle below southern Spain. Synthetic slab tests In order to tackle the resolution problem from a different angle let us consider the second question: if a subducted slab exists in the upper mantle below the Betic and Alboran Sea, can it be imaged with sufficient reliability using our data set? We investigate this with sensitivity analysis using synthetic slab models. We constructed a synthetic velocity anomaly model of a subducting slab that resembles the

Ll’fHOSI’HERt

Fig. 12. Average

station-delays

than 20 observations. the symbols plotted. the station-delay

BELOW

27

S SPAIN

for regional

The stations

stations

are located

The size of the symbols

(see key). The station-delays

from delay times relative

with more

at the center is proportional

of to

are computed

to model PM?.

shape and amplitude (2% in the core and tapered to 0% at the edge) of the anomaly mapped in IPl. The spatial outlines of the synthetic model are displayed in Figures 13 and 14. Synthetic delay-times are computed from this model using the actual ray paths. We performed an inversion of the exact synthetic data (result SYNO, data fit 60%), an inversion of noisy synthetic data (result SYNI, signal to noise ratio 0.15, data fit 4%), and an inversion of noisy synthetic data with epicentral distances between 0 and 30” (result SYN2, data fit 4%). in map-view we only show result SYNl (Fig. 13). From SYNl we can readily infer that, even in the case of a very poor data fit, the synthetic model is well recovered. Al1 anomaly patterns outside the boundaries of the synthetic slab are of course resolution artifacts. Specificaliy, anomaly leaking into the lower mantle is clearly visible. Cross-sections through the sensitivity test results are displayed in Figure 14. Note that even in the case SYNO (no data noise!) small resolution problems exist. The smearing of the synthetic slab anomaly along the main directions of ray-illumination corroborates our earlier findings in

M.J. BLANC0

1Synthetic

slab: SYNI

Fig. 13. Tomographic

result SYNI of a noisy data inversion

shape

and amplitude

indicated -0.5

the slab-like

with a thick white

and 1.5V. Anomalies

anomaly

1

imaged

line visible in the panels outside

line are all imaging

+ 1.5%

-0.5% for a synthetic

beneath

slab model. The synthetic

the Betic and Alboran

with depths artifacts.

AND W. SPAKMAN

between

area.

247 and 635 km. Notice

Note that the core of the synthetic

SYNI to IP1.

slab is designed

The location

to resemble

of the synthetic

in

slab is

that the contouring

limits are

slab is well recovered.

Compare

P-WAVE

VELOC‘ITY

STRUCTURE

OF THE MANTLE:

EVIDENCE

FOR SURDUCTED

I.1 I‘HOSPHEKt

BtLOW

29

S SPAIN

200 400 g

600

f &? 800 0 1000 1200

200 400 $

600

8

800 1000 1200

200 400 g

600

$

600 1000 1200

-0.5%

+I 25%

Fig. 14. Two cross-sections results.

The thick white

through line indicates

recovered

in SYNO (from

epicentral

distances

noise-free

the synthetic the actual data)

less than 30”). Contouring

and

results

SYNO, 1, and 2. The cross-sections

boundary SYNl

of the synthetic

(from

as in Figure

SYNO and SYNl should be compared

+I .5%

-0.5%

are the same as for the actual

slab. Note that the core of the synthetic

very noisy data)

and

13. Notice that imaging

poorly

artifacts

in SYN2 resemble

to result IPI. Result SYNZ must be compared

(from

data

slab is well

very noisy data

with

those in the real data results. to IP3.

30

M.J. BLANC0

AND W. SPAKMAN

200 400 600 f $ 0

800 1000 1200

200 400 z f B

600 800 1000 1200

200 400 g

600 800 1000 1200

Fig. 15. Two cross-sections through the synthetic results SYN4, SYNS, and SYN6 along the same profiles as in Figure 14. The boundary of the synthetic (input) model is indicated by thick white contour lines and consists for SYN4 of a high velocity lithosphere (+ 2%) of 70 km, a low velocity zone ( - 2%) between 70 and 220 km, and below 220 km the model is the same as in the first slab test (SYNl). For SYN5 the lithosphere is 120 km in thickness. SYN6 is made from a continuous slab between 0 and 670 km and - 1% low velocities surrounding the slab between 0 and 220 km. Contouring as in Figure 13. Compare SYN4, SYNS, and SYN6 to SYNI and IPI.

P-WAVE

VELOCITY

STRUCTURE

OF THE MANTLE:

EVIDENCE

FOR SUBDUCTED

the data results IPl, IP2, and IP3. From the result SYN2 we can infer that short-distance data primarily reconstruct the upper part of the synthetic slab anomaly. The occurrence of strong lateral smearing in result IP3 is well supported by the similar resolution artifacts in SYN2. In spite of the poor resolution for small objects (cf. spike test, Fig. 51, the tests demonstrate that the resolution provided by the reference-model rays in mapping larger objects is reasonably good beneath southern Spain, even when we force a small signal to noise ratio on the synthetic data. Therefore, if a slab-like structure is present below southern Spain and the Betic, we can detect it with sufficient reliability. Moreover, from the close correlation between the anomaly smearing in the synthetics (SYNl, SYN2) and data results (IPI, IP3) (where the “slab anomaly” is concerned) we find another (but indirect) indication that a slab-like velocity structure actually exists at depths between 200 and 700 km below southern Spain and the Alboran sea. The real data results IPl and IP3 suggest that the slab structure is overlain by low velocities and hence may not be connected to the Betic-Alboran lithosphere. Similar observations have been made for the Hellenic slab below western Greece and for subduction beneath Italy; they are interpreted as an image of slab detachment from the shallow lithosphere (Spakman et al., 1988; Spakman, 1990). To obtain insight into the resolution for a slab detachment configuration, we performed additional tests with three synthetic slab models, modelled below the Betic-Alboran region. The first model possesses a high-velocity lithosphere of 70 km, a low-velocity zone of - 2% from 70 to 220 km, and, for depths larger than 220 km, the model is the same as before. The inversion result is called SYN4. The second model has a thicker lithosphere of 120 km (+ 2%), and a thinner low-velocity zone from 120 to 220 km (-2%). This inversion result is called SYNS. The last model consists of a continuous slab from 0 to 670 km of a +2% anomaly value, which is surrounded by a low-velocity lithosphere and asthenosphere of - 1% from 0 to 220 km, for the entire Iberian region. The inversion result is called SYN6. In all these tests we added normally dis-

LITHOSPHERE

BELOW

S SPAIN

31

tributed noise with a standard deviation of 0.5 s to the synthetic data. The three results are displayed in Figure 15 for the two cross-sections and can be compared to IPl, and SYNl. We can infer that the slab detachment configuration is resolvable (SYN4), even for a thin low-velocity zone (SYNS), although in the latter case some smearing across the low-velocity zone does occur. The test SYN6 is made to demonstrate that we could image a continuous slab, even if the upper 220 km of the slab is surrounded by only negative anomalies. This demonstrates that low velocities above the slab image in IPl (i.e., above 200 km) are probably due to other causes than anomaly smearing from adjacent low-velocity structures (if these were present). Below 200 km the anomaly patterns in SYN4, SYNS, and SYN6 are quite similar to those of SYNl (Fig. 14) which indicates that the large scale structure of the lithosphereasthenosphere system in the Betic-Alboran region is to a large degree independently resolvable from the lower part of the upper mantle. Since the slab image in IPl resembles more the detachment configuration (or a configuration where a high-velocity lithosphere seems largely absent locally), we use these tests.in support of our interpretation of slab detachment below the Betic-Alboran region. Discussion

and conclusions

In this paper we applied delay-time tomography to retrieve the P-wave velocity structure of the mantle beneath the Iberian Peninsula. The method of tomography applied attempts to account for the effects of event mislocation, origin time errors, and station statics on the data. The data are well centered with respect to the traveltime predictions determined from radial velocity model PM2 (Fig. 3) which at least indicates that the non-linear effects of reference model problems related to ISC delay-times (Van der Hilst and Spakman, 1989) have been reduced. Other effects of non-linearity (e.g., ray bending due lo the velocity heterogeneity) are not considered in the present method, nor the effects of possibly remaining systematic errors/signal in the data (see also below). Therefore, the results obtained

32

here can only be considered as a first approximation to the real Earth structure, although we do not anticipate a significant improvement of the results when we would account for ray bending effects in a more elaborate approach to delay-time tomography. This, because the noise level in the data is rather high and the resolving power of the ray set suffers from dominating directions in ray illumination as inferred from the hit-count patterns. An improvement could be obtained by including data from new station-event combinations allowing for a more isotropic illumination of the mantle region below the Iberian Peninsuia. For other uncertainties and problems in delaytime tomography see Spakman (1991). All our inferences about resolution are derived from sensitivity analyses with synthetic velocity models. It is important to realize that sensitivity analysis gives only insight in the resolving power of the ray set employed, and not of the actual delay-time data and real ray paths. Care is taken to mimic with the sensitivity tests the convergence behaviour of the real data inversions (by adding normalty distributed noise with a sufficiently large standard deviation to the synthetic data). Sensitivity tests are very useful in the case that the reference ray paths are good approximations of the actual ray paths, and if the noise in the actual delay-times is uncorrelated, has zero mean and possesses a symmetric distribution. As remarked above, we believe that the only ray path deviations of the real paths with respect to the reference paths are primarily due to the velocity heterogeneity since we corrected for reference model influences. The large scale velocity heterogenei~ itself is of modest amplitude (a few percent at most) and will therefore not result in too large ray path deviations. Our ray set is probably accurate enough to infer resolution estimates from sensitivity analyses in this respect. With regard to the noise in the delay-time data, we can say that the data between -3 and +3 s are about symmetrically distributed and therefore it is very reasonable to assume that also the data errors will be symmetricahy distributed. Systematic and correlated errors errors remain difficult to assess. One test is useful to investigate effects of noise with the same statistical distribu-

M.J. BLANC0

AND

W. SPAKMAN

tion as the data: the permuted-data test (for details, see Spakman and Nolet, 1988; Spakman, 1991). In this test, the actual delay-time data are randomly permuted prior to inversion while keeping the order of the tomographic equations fixed. In this way the correlation between delay-times on the one hand and ray paths, events and stations on the other hand is destroyed, but the statistical distribution of the delay-times stays the same. The inversion of the permuted data set rendered only a few small scale anomaly patterns with small amplitudes on the order of O.l-0.3% and with no correlation with the patterns in result IPl, which demonstrates that large errors (with data’s distribution) will not lead to large scale and/or high amplitude heterogeneity in the inversion, an observation which gives additional credit to the applicability of sensitivity tests for resolution analysis. Furthermore, the result of this test leads us to conclude that IPl is a significant inversion result on the structure of the Iberian mantle (see Spakman, 1991, for how to arrive at this conclusion). In fact, this is the only indication we have that IPl is a significant result. Sensitivity tests are not informative in this respect. From the start we have concluded that the resolving power of data set for imaging small-scale structure on the order of 100 km is very limited, except for parts of the upper 150 km. From the block-model inversion (Figs. 6, 7) we inferred reasonable resolution for mantle structures with dimensions of a few hundred km or more. In view of the poor resolution for detail our analysis has been focussed on locating the spatially large velocity structure that has caused the mapping of a slab-like anomaly in the upper mantle beneath the eastern Betic and the Alboran sea. The block inversion result already indicated the good potential of the data set to resotve the slab structure. We demonstrated from the results of inversions of subsets of the data for different epicentral distance ranges that the likely location for this velocity structure is at least in the upper mantle below southern Spain. A lower mantle extension of the slab anomaly is not excluded by this approach. To arrive at this point we also invoked the argument that the event mislocation and sta-

P-WAVE

VEL.OCITY

STRUCTURE

OF THE

MANTLE:

EVIDENCE

FOR

tion parameters largely absorb the effects of lateral heterogeneity outside the cell model. We proceeded our analysis with a resolution test for a synthetic slab structure located at the same place and with a similar amplitude as the slab-like anomaly in the data inversion results. Our conclusion from this test is that the resolution for large objects below southern Spain is reasonably good. We put effort in demonstrating how image blurring effects are related to the hit-count pattern. Specifically, the similarity between the resolution artifacts (e.g., anomaly smearing) in the real data results and in the synthetic results point (indirectly> at the actual existence of a high-velocity body below southern Spain. In conclusion and adding all inferences, we have demonstrated that it is likely that a slab-like structure resembling the geometry as inferred from the inversion rest&s actually exists beneath southern Spain. We interpret the anomaly between 200 km and 700 km as the mapping of a subducted slab. From the slab image in IPl (Fig. 10) we cannot infer a clear dip-direction more accurate than about vertical. In map-view the slab has a clear SW-NE strike direction and is of limited lateral extent. The slab is located primarily below the Alboran sea and eastern Betic. The slab-inte~retation is corroborated by the existence of very deep seismic&y at 640 km located in the core of the slab anomaly. An indication for a subducted slab was found earlier by Spakman (1986, 1990) who obtained a slab image from inverting primarily short-distance data. The image is poorer resolved but is similar in shape to that in result IP2. We are not the first to suggest that subducted lithosphere is present below southern Spain. Its presence has been invoked by many authors to explain the deep seismicity (e.g., Udias et al., 1976; Chung and Kanamori, 1976; Buforn et al., 1988, 1991). Our reasoning, however, is opposite: we have found slab of considerabIe length and merely invoke the presence of the deep seismicity to corroborate this result. Our results suggest that the slab is not connected to high-velocity lithosphere at the surface and can be interpreted as detached. The resolution provided by our ray set has been demon-

SUBDUCTED

LITHOSPHERE

BELOW

S SPAIN

33

strated to be sufficient for imaging slab detachment and our results resemble more a detached slab configuration than that of a continuous slab. In fact, it seems that locally also a distinct highvelocity lithosphere is lacking in the Betic-Alboran region, a result which has already been inferred from surface wave studies (Panza et al., 1980; Marillier and Mueller, 1985; Suhadolc and Panza, 1989). The low velocities found above the slab may be partly due to continental lithosphere involved in the present interaction between Iberia and the African plate. Another cause of the lowvelocity material found directly above the slab may be asthenospheric material which by lateral inflow has occupied the space created by detachment. At the time it occurs an important consequence of slab detachment is the dynamic response of the crust-lithosphere system directly above the slab. Prior to detachment the elastic energy stored in the down-bended part of the slab above the future detachment level will be released upon slab detachment and may cause large uplift in the region above the slab (Wortel and Spakman, 1992). In the Betic region the most recent phase of important uplift occurred in the early Miocene. To explain this tectonic phase Platt and Vissers (1989) postulated the detachment of a “cold lithosphere root” below the Betic region. Combining our inferences with those of Wortel and Spakman, and of PIatt and Vissers we can tentatively place the time of slab detachment in the early Miocene and consequently subduction took place during at least part of the Oligocene,

This research was carried out while M.J. Blanc0 was a visiting scientist at the Department of Geophysics of the University of Utrecht (fan 1989) supported by a grant from the Universidad Complutense (Madrid). We have benefited from the constructive comments made by E. Kissling, M.J.R. Wortel and an anonymous reviewer. References Badal, J., Corchete, V., Payo, G., Canas, J.A., Pujades, L. and Ser6n, F.J., 1990. Processing and inversion of long-period

M.J. BLANC0

34

surface-wave

data collected

in the Iberian

Peninsula.

Geo-

phys. J. Int., 100: 193-202. Banda,

E., Ansorge,

Structure

of the

Balearic

Islands

J., Boloix,

M. and

crust

upper

and

(western

Cordoba, mantle

D., 1980.

beneath

Mediterranean).

Earth

the

Planet.

E., Suriiiach,

la Parte, sot.,

E., Aparicio,

A., Sierra,

E., 1981. Crust and upper

central

Iberian

Meseta

J. and Ruiz de

mantle

(Spain).

structure

Geophys.

of the

J.R.

E., Udias,

Gallart,

A., Mueller,

Spain

M.,

structure

be-

deep

seismic

experiments.

Inter.,

31: 277-280.

from

Phys. Earth Buforn,

J., Boloix,

A., 1983. Crustal

J. and Aparicio,

neath

Planet.

E., Udias,

St., Mezcua, sounding

A. and Mezcua,

focal mechanisms

J., 1988. Seismicity

in south Spain.

Bull. Seismol.

Sot. Am.,

i

E., Udias,

deep

A., Mezcua,

earthquake W.-Y.

tectonic quake

under

J. and Madariaga,

south

Spain.

R., 1991. A

Bull. Seismol.

implications

H., 1976. Source

of the

Spanish

process

deep-focus

29, 1954. Phys. Earth

Planet.

M. and

structures

13:

Trampert,

J., 1989.

Large-scale

in the Euro-Mediterranean

P-velocity

area.

Geophys.

J.

E., 1988. Geotomography

Rev. Geophys., Marillier,

F. and

ranean

Mueller,

region

plates.

systems.

J. Comput.

structure

recording Paige,

under

Trans.

linear equations

Math.

rope Appl. Paulssen, Iberia

Geophys.,

of autonomously

1991. Delay

below

nor. Geophys. Spakman,

and sparse

least squares.

surface

waves

inferred

from P-wave

coda.

(Editors),

Seismic

Tectonophysics,

Pure

means

of a long period Sot, 20: 493-508.

J. and Vissers,

array.

Geophys.

by J.R.

in the Iberian

of the seismicity

Penin-

in this area.

Sot, 30: 85-99. 1989. Extensional

lithosphere:

Sea and Gibraltar

collapse

arc. Geology,

of

hypothesis 17: 540-

and

and

M.J.R.

and

zone:

implications.

Miaccu-

In: N.J.

S.A.L.

Cloetingh

a Survey

of Recent

Geodynamics.

W., Stein,

Vlaar,

S., Van

Zone Tomography.

N.J.,

a tomographic Geophys.

duction

Reidel,

der

1988. The

image

Res. Lett.,

Hilst,

R.D.

experiments

and

and

Wortel,

for NW Pacific

Geophys.

its

15: 60-63. Sub-

Res. Lett., 16: 1097-

1100. W., Van der Lee, S. and Van der Hilst, R.D., 1991.

Structure depth

of the European-Mediterranean

of 1400 km: preliminary

Suhadolc,

Terra

Abstr.,

P. and Panza,

physical

mantle

results

to a

from P delay-time

3(l): 173.

G.F., 1989. Physical system

data. In: Boriani

properties

in Europe

et al. (Editors),

Naz. Lincei,

A., Lopez-Arroyo,

Rome,

of the

from

geo-

The Lithosphere

pp. 15-40.

A. and Mezcua,

of the Azores-Alboran

J., 1976. Seismo-

region.

Tectonophysics,

31: 259-289. R.D.,

Delay-Time Structure

1990. Tomography

Data

and

below

Utrecht,

Utrecht, R.D.

the reference

the

Res. Lett.,

and Spakman,

structure

Modeling

tomographic

Univ. of

inversions:

plate. Geophys.

W., 1992. Structure

lithosphere

Proc. K. Ned. Akad.

and

Mantle

Thesis,

W., 1989. Importance

below the Caribbean

A., Spakman,

velocity

region.

in linearized

in the

imaging

beneath

G. Nolet

G., 1989. A reference

of the upper

Europe.

(Editors),

of the Lithosphere.

and dy-

Mediterranean

Wet., 1992: in press.

W. and Nolet,

model for tomographic

340.

Caribbean

and Spakman,

model

of subducted

region,

with P, PP and pP

Three-Dimensional

16: 1093-1097.

M.J.R.,

namics

the

250 pp.

images of subduction

Panza

A working

upper

algorithms,

tomography.

Wortel

in Seismology

1989. Resolution

Wortel, mantle

time

Geophysics:

M.J.R.,

Zielhuis, velocities

R.L.M.,

continental

J. Gallart

of the Iberian

and upper

triangular

implications

J.R. Astron.

Studies

in

221: 111-123.

Payo, G., 1972. Crust-mantle sula and tectonic

structure

In: J. Badal,

of the crust

Astron.

for the Alboran

waves.

crustal

M.J.R.

subduction

Van der Hilst,

J., 1993. The

G., 1970. Structure

thickened

body

in Eu-

118: 1209-1213. Visser,

Geophys.

and

G., 1988. Imaging

W., Wortel,

Van der Hilst, system

of the

and Asia

pp. 155-188.

tectonics

G., 1980. The gross

of central

Nova, 2: 542-533.

tomography

in delay

Mathematical

Hellenic Spakman,

mantle

Terra

the Mediterranean,

resolution

G. Nolet,

Spakman,

53, 200

J. Int., 107: 309-332.

in Italy. Acad.

ACM

65: 145-

Delay Time Tomography.

time

Europe,

W. and Nolet, and

Vlaar,

Udias,

an algorithm

Res

in connec-

Geol. Ultraiectina,

of the upper

lithosphere-asthenosphere

and two-di-

1982. LSQR:

H. and

Peninsula.

noisy

Res., 95: 8499-8512.

St. and Calgagnile,

seismic

and H. Paulssen

543.

inversion

of the lithosphere-asthenosphere

from

and

Soft., 8: 43-71

G.F., Mueller,

features

Platt,

W.,

tomography.

118: 113-130.

the network

M.A.,

zone between

Phys., 61: 463-482.

J. Geophys.

C.C. and Saunders,

Panza,

data.

Mediter-

inadequate

waveform

seismographs.

for sparse

Payo,

or resolving

G., 1990. Partitioned

mensional

transition

Tectonophysics,

1985. Solving

tomographic

St., 1985. The western

Mantle

and the Mediterranean.

Spakman,

as an upper-mantle

two lithospheric G.,

with local earthquake

26(4): 659-698.

Eurasia

Geol. Mijnbouw,

Utrecht.

W., 1990. Images

Europe

Dordrecht,

Int., 99: 583-594.

Nolet,

PP. Spakman,

and earth-

J. Geophys.

beneath

Tethys.

W., 1988. Upper

geodynamic

Granet,

Nolet,

Spakman,

(Editors),

Inter.,

to surface

153.

Sot.

85-96.

I&sling,

W., 1986. Subduction

Developments

and Kanamori,

of March

Spakman,

racy

Am., 81: 1403-1407. Chung,

and the Mediterranean.

of surface

II: Application

93: 12067-12080.

Spakman,

and

inversions

waves in Europe

mantle

78: 2008-2024. Buforn,

scale waveform

heterogeneity,

Thesis, Univ. Utrecht,

Astron.

67: 779-789.

Banda,

R., 1988. Large

waves for lateral

tion with the Mesozoic

Sci. Lett., 49: 1219-230. Banda,

Snieder,

AND W. SPAKMAN

Digital Plenum,

mantle

In: R. Cassinis, Seismology London,

shear G. and

pp. 333-