3. Tomography Study
The data set used for this study was selected from data recorded by multiple networks in Spain and Morocco (Figure 3). Most of the 96 regional stations from which data were used are predominantly short-period vertical component with a limited number of three-component and intermediate band stations. Digital waveforms from approximately 200 local/regional and teleseismic events between 1989 and 1996 were provided by Instituto Geografico Nacional (IGN) in Madrid (Spain); Instituto Andaluz de Geofisica (IAG) in Granada (Spain) and Centre National de Coordination et de Planification de la Recherche Scientifique et Technique (CNCPRST) in Rabat (Morocco). Supplementary analog records were also provided by CNCPRST and Physique du Globe at Mohammed V University (MOH V) in Rabat. P and S arrivals were repicked for the local events and P arrivals from the teleseismic events (i.e., between 30 and 90 degrees epicentral distance). These picks were associated and merged with catalog data from the International Seismological Centre (ISC) (1964-1994), USGS National Earthquake Information Center (NEIC) Earthquake Data Report (1995-August 1998), and Moroccan CNCPRST bulletins to generate comprehensive databases of local and teleseismic arrivals recorded in Morocco, Spain, Portugal, and Algeria from 1964-August 1998. Errors in the central clock for CNCPRST network required that readings from this network not be merged with other networks during certain periods. The combined local/regional data set contained 7100 events with 66,200 P and 48,600 S arrivals. The S arrivals for local/regional events were included to improve constraints on the earthquake locations (particularly depth). The teleseismic data set comprised 15,700 events with 47,700 P arrivals. Teleseismic arrivals reported as emergent were not extracted from the ISC and USGS catalogs, as these are likely to contain significant uncertainties.
The independent repicking of phases by the authors using the waveform data allowed some quantitative
assessment of the reading uncertainty (sreading)
and hence determination of reasonable data weights (1/s2reading) for the larger catalog data set. Differences between the independent picks (assuming
no systematic biases) made by the authors and those reported in bulletins were calculated, and the standard deviation
of this difference (SDD) was determined for the different reported qualitative picking precisions (e.g., impulsive
as opposed to emergent) for each network. The SDD of measurements was calculated after removal of any arrival time
differences (outliers) that exceeded a certain threshold value. A threshold of 3 s was found to significantly reduce
the SDD. This threshold was included to simulate the effect of strongly down weighting travel time residuals greater
than a user-defined cutoff in the velocity inversion (discussed below). A 3-s residual cutoff still allows residuals
resulting from significant differences between the real Earth and the model (e.g., an average 3% higher velocity
relative to 8.0 km/s over 800 km) to contribute to the inversion. The standard deviations calculated after a 3-s
cutoff are summarized in Table 1. Not surprisingly, arrivals
reported as impulsive in the ISC and USGS bulletins are of better quality than those reported as emergent especially
for the Spanish networks where picks are made from digital records. Teleseismic residuals were also examined and
found to be consistent with clear azimuthal variations for most stations in Spain. However, some stations in Morocco
showed some scatter and were assigned a larger reading uncertainty. An additional uncertainty of 0.1 s was also
added to account for inaccuracies resulting from the need to approximate the Earth structure using discrete blocks
of constant velocity and inexact ray tracing. Local and teleseismic events recorded by six or more stations in
the Alboran region were extracted from the database resulting in a final selected data set of 4500 local events
with 55,800 P arrivals and 39,500 S arrivals (Figure 4a) and 1600 teleseismic events with 16,400 arrivals (Figure 5).
Seismic source regions that provide local data are scattered, indicating the distributed nature
of active deformation in the region (Figure 4). However, the
majority of well-constrained crustal earthquakes are concentrated in the Betic and Rif mountain belts. Intermediate-depth
earthquakes occur predominantly in two regions : (1) A diffuse zone off the coast of Spain in the Gulf of Cadiz
and (2) relatively concentrated N-S line spanning the Alboran Sea and dipping south from crustal depths beneath
the Betics to a depth of ~150 km beneath the center of the Alboran Sea. This dipping band of seismicity is overlain
by a relatively aseismic zone. A limited number of well-constrained intermediate-depth earthquakes have also been
located off the Atlantic coast of northern Morocco . No well-constrained earthquakes have been located between
depths of 200 and 550 km in the region. However, four events have been located at depths of 550-650 km beneath
southern Spain . One of these events (not shown in Figure 4),
which occurred in 1954 with Mw = 7.8 , is the largest instrumentally recorded earthquake in the region.
The teleseismic events are located predominantly along the subduction zones of the Americas, the Atlantic ridge system, Zagros and the India-Asia collision zone (Figure 5). Azimuthal coverage is generally good, especially from the east and west, although there is limited illumination from the south.
3.2. Relocation of Local and Regional Events
An initial one-dimensional (1-D) velocity model (Figure 6 and Table 2) was constructed from a reference model (PM2) determined for Europe by Spakman et al. and regional refraction surveys . Each event was relocated using the same starting epicentral location and origin time but multiple starting hypocentral depths of 0.5, 14, 40, 90, 300, and 600 km. The depth was fixed for the early location iterations to reduce the influence of the starting epicentral location and origin time on the final hypocenter location. Residual cutoffs (used to downweight outliers) were progressively reduced from 60 to 3 s. Despite the different starting depths, the solutions for a single event usually converged to a single hypocentral location and origin time. However, approximately 10% of events converged to two or sometimes even three locations, suggesting the presence of multiple residual minima. In these cases, the location with the lowest weighted RMS residual was chosen. Events with poorly constrained locations were removed using criteria 1-8 shown in Table 3. To ensure that only events with locations that were relatively insensitive to the assumed Earth model were used for the inversion, the events were relocated in 10 different velocity models generated by adding a random perturbation (uniformly distributed between +/- 5%) to each layer of the reference 1-D model. Any event that moved by greater than 15 km horizontally, 20 km vertically, or 25 km in total was discarded . Finally, to improve homogeneity of the ray coverage in the model and reduce the biasing of station corrections and 1-D inversions by arrivals from aftershock sequences, the local data set was geometrically thinned by preserving the event with the most recorded arrivals in each 5 km cube of the model . A final local data set of 1100 events consisting of 16,000 P and 11,200 S arrivals was used for subsequent 1-D inversions. The three deep earthquakes with available data were not included in this data set owing to excessive trade-off between event depth and the velocity model (greater than 20 km variations).
3.3. Inversion Method
The 1-D and 3-D inversions were conducted using the inversion code SPHYPIT90 developed by Roecker and further modified by numerous researchers. The algorithm follows the nonlinear maximum likelihood formalism of Tarantola and Vallette to iterate toward the bestfit solution for the model parameters:
pk+1 = pk+(GkTCdd-1Gk+Cpp-1)-1[GkTCdd-1Gk[d-g(pk)]-Cpp-1(pk-p0)], (1)
p0	a priori or starting estimate of the model parameters;
pk	estimate of model parameters for the kth iteration;
d	data vector of observed travel times;
g	expected travel times (from ray tracing the model);
Gk	matrix of partial derivatives (dg/dp);
Cdd	diagonal a priori data covariance matrix of reading uncertainties;
Cpp	a priori model covariance matrix.
The elements of Cdd are modified by two dynamic weights scaled between 0 and 1. Local/regional events are multiplied by a distance weight to account for a decrease in picking accuracy at longer distances owing to a reduction in the amplitude and frequency content of phases with distance. Local and teleseismic readings with relative residuals greater than 3 s are also significantly down-weighted to reduce the effect of outliers. So the elements of Cdd are
Cdd(i) = s2reading/(distance_weight x residual_weight).	(2)
Cpp is our estimated uncertainty in the slownesses assigned to the a priori
model parameters p0 (1-D inversion model). In this case, Cpp is a diagonal
matrix with no covariance assumed between model parameters. The elements of Cpp were chosen to
correspond to a 0.2 km/s uncertainty in the a priori model. Other values were investigated, including allowing
larger uncertainties in the crustal layers; however, a constant uncertainty of 0.2 km/s was found to be optimal
with a minimum of unconstrained assumptions.
Travel times and partial derivatives are recalculated for each iteration by ray tracing through an average 1-D model calculated between each station-event pair and accumulating the travel time along this path in the 3-D velocity model. This approach is not true 3-D ray tracing, but it was found to be a good approximation by Thurber and Ellsworth . Teleseismic arrivals are included using a modification of the Aki, Christoffersson, and Husebye (ACH) method by ray tracing to the base of the model volume using the radial IASPEI91 model before tracing onto the station through the 3-D model. Parameter separation allows the new slowness model and event relocations to be determined in separate steps. First, the slowness inversion is conducted after the removal of the hypocenter effects, then the local/regional events are relocated in the new velocity model using a separate location program. After each slowness-relocation iteration, events that fail criteria 1-10 shown in Table 3 are removed. The matrix inversions were performed using singular value decomposition. Although we inverted for both P and S wave velocity structure in the lithospheric layers, we will only discuss the P inversion results below. As mentioned previously, the main reason for including local/regional S wave arrivals was to provide additional stability to the earthquake locations.
3.4. One-Dimensional Inversion
The local data set was used to invert for a bestfit 1-D model for the upper 150 km. Station corrections were calculated before each iteration step from the mean travel time residual at a station. The solution converged after three iterations. The bestfit 1-D model is shown in Table 2 and in conjunction with the station corrections results in a residual variance improvement of 30%. Mean station delays show a strong regional variation (Figure 7), suggesting that some crustal structure is being mapped into the station corrections. The most significant features are the consistently late arrivals to the stations in the Strait of Gibraltar region and the variation in the stations east of Malaga from early in the south to marginally late in the north. The substantial station delays in the Gibraltar region may be partly caused by the significant thickness (up to 10 km) of intensely deformed low-velocity (2.3-5.2 km/s) sediments that were observed in the flysch regions by refraction experiments . The significant shortening and associated faulting that occurred in this region probably further contributes to the low velocities in this area. The high-velocity signal in the station corrections on the southern coast of Spain was also identified by refraction studies and is discussed in more detail below. These station corrections were held constant for the 3-D inversions.
3.5. Three-Dimensional Inversion
The 3-D model is parameterized by subdividing each layer into blocks using vertical interfaces. After a number of trial inversions and examination of the ray coverage, we chose the final model parameterization shown in Figure 8. The upper four layers were divided into 50-km blocks to take advantage of the dense ray coverage from the local/regional events. The deeper layers that are controlled predominantly by teleseismic data were divided into 100-km blocks. Figure 9 shows a hit map for a 50-km block parameterization throughout the model volume, to provide more detail of illumination in the deeper layers.
Four iterations were performed with the last iteration providing an improvement in fit of less
than 1%. Comparison of the images after the first and fourth iterations indicated that first-order structures are
visible after only one iteration but significant changes are obtained with further iterations, suggesting that
the inversion is not completely linear. As a metric, we compared the diagonal elements of the partial derivative
matrix Gk for the first and last iterations. These elements should remain constant in a linear
inversion. Although the majority of elements changed by less than 3%, some fluctuated in value by up to 10%, further
suggesting that it was prudent to use a nonlinear approach. The final 3-D velocity model produced an overall variance
reduction of 32% relative to the 1-D model. The final variance was 0.35 s2, whereas the a priori variance
was 0.19 s2, suggesting that the model was not overparameterized.
As the location of block boundaries is arbitrary and ray illumination is not homogeneous, sharp velocity contrasts that do not lie close to a block boundary may be smoothed over a block or displaced ("leaked") to a block boundary. Furthermore, any artifacts that are introduced by a particular choice of block boundaries will probably not be detected by the standard tests used to assess the quality and stability of models. In an attempt to limit these problems, we used a block shift-and-average approach . One hundred separate inversions were conducted using the same starting 1-D model but with the block boundaries uniformly shifted by (10n+5) km in the north-south direction and (10m+5) km in the east-west direction for n = -5 to 4 and m = -5 to 4. The mean RMS residual for all these different model parameterizations was 0.349 s2 with a 2s of 0.006 s2 (minimum = 0.345 s2 and maximum = 0.368 s2 ), suggesting that no one model parameterization provide a better (in a statistically significant way) fit than any other. All 100 models were resampled onto the same grid with a 10-km horizontal node spacing. The mean velocity at each grid node was calculated resulting in a smooth model that was not a function of a particular parameterization with any significant anomalies in the final model caused by an anomaly being imaged by many different parameterizations. A model derived using only 16 block shifts of 25 km showed identical features (allowing for sampling differences) so this less computationally demanding shifting scheme was used for the synthetic tests described below. Some measure of the parameterization dependent stability of the velocity model was provided by the standard deviation of the velocities. Areas of high variance will not only occur in regions of poor stability but may also indicate the presence of large but spatially limited (on the scale of blocks) velocity anomalies. An alternative approach for limiting the effect of an arbitrary parameterization, might be to use variable block sizes with smaller blocks in regions of good illumination .
3.6. Uncertainty Analysis
The influence of systematic and random errors on the model is hard to quantify owing to the nonlinear nature of seismic tomography, the difficulty of measuring the amount of noise in the input data, and the effects of parameterization. A number of traditional and nontraditional methods of uncertainty estimation were used to ensure that only robust well-constrained features were interpreted. These methods are reviewed below using an example cross section located 25 km north and parallel to cross section BB' (Figure 10). Views of the entire model volume for the different tests summarized in Figure 10 and discussed below may be found on our website at atlas.geo.cornell.edu (You are here). Figure 10a shows a cross section through the final shift-and-average velocity model generated by stacking the results of 100 inversions using different model parameterizations. Figure 10b shows the velocity model resulting from a single inversion demonstrating that the first-order features are the same as those found in the shift-and-average model.
3.6.1. Bootstrap resampling of residuals. The increased speed of modern workstations now
allows the use of nonlinear uncertainty estimators to evaluate inversion results. Hearn and Ni used a bootstrap
resampling technique to determine uncertainties for 2-D Pn velocity models. The same approach was used to
evaluate the 3-D tomography model determined in this study. The bootstrap method involves repeated inversions of
resampled data sets and examining the variance in the derived models. To generate a resampled data set, arrivals
times are randomly sampled from the data pool. However, after a data point is sampled it is returned to the pool
so it may be picked multiple times. The same number of data are picked from the pool as existed in the original
data set but each data set will be different owing to the aforementioned data redundancy and resulting omission
of other data points. A velocity model is determined from these resampled data using the same inversion parameterization.
This model is stored and the process repeated approximately 100 times to generate robust parameter distributions.
The standard deviation of a particular model parameter is a measure of the uncertainty in that value in the final
model. An important strength of a bootstrap method is that, in contrast to Monte Carlo simulations, it does not
require an estimate of the noise in the data to assess the uncertainty in the resulting model. The resampling technique
allows the noise in the data to contribute to the variance in the derived models. However, it should be noted that
in this case the bootstrap approach does not remove the need for noise estimates to measure model uncertainty because
they are required by the inversion method (the elements of Cdd). A bootstrap conducted by resampling
events rather than individual phases produced virtually identical results. The bootstrap uncertainties are invariably
higher than standard errors calculated using the a posteriori covariance matrix (compare Figures
10c and 10d). The standard errors are less than 1% beneath the Betics in the lithospheric layers
and less than 1% beneath the Alboran, Betics, and Rif in the deeper layers. The bootstrap uncertainties determined
for layer depths of 100-670 km, which are controlled predominantly by teleseismic data, were encouraging with important
features imaged with anomalies greater than twice the bootstrap error (see Figures
10b and 10d). Although the crust and uppermost mantle layers showed significant variance (sometimes
approaching 4%), these regions also contained the largest perturbations with robust anomalies again above twice
the bootstrap uncertainty. A useful method for interpreting our model data was to subtract between 1 and 2 times
the bootstrap uncertainty from the imaged anomalies setting an anomaly to zero if the anomaly changed sign after
this subtraction. Only the more robust anomalies remained for interpretation.
3.6.2. Resolution matrix. The resolution matrix provides an estimate of the relative contributions
of the observations and the a priori model to the final model. The diagonal elements of the resolution matrix are
very close to one beneath the Betics, suggesting that the ray geometry and choice of model parameterization and
damping parameters allows the imaged crust and mantle structure to be controlled by the observations rather than
the a priori model (Figure 10e). Crustal structure beneath
Morocco is not well controlled owing to the limited number of local earthquakes. Off-diagonal elements of the resolution
matrix were also examined to evaluate smearing. The off-diagonal elements are usually an order of magnitude less
than the diagonal elements suggesting that at least in theory there is very limited smearing of anomalies or generation
of compensating anomalies in adjacent blocks.
3.6.3. Spike test. Spike tests are traditionally used to assess resolution in tomographic
inversions. Synthetic data are generated by ray tracing a velocity model consisting of a regular 3-D pattern of
perturbations to the reference model (e.g., Figure 10f).
A representative amount of noise is often added and the synthetic travel times inverted following the same procedure
that was used for the real data. The recovered velocity model is compared to the velocity model used to generate
the synthetics. Regions where the recovered model closely matches the input model are considered well resolved
in the real velocity model, whereas a poorly recovered region is considered to be suspect in the final model. The
addition of noise usually degrades the recovered model, but in this case it has a limited effect (compare Figures 10g and 10h).
The degree of recovery is also sensitive to the geometry of the synthetic test and in particular its relationship to the model parameterization used to conduct the inversion. Figure 10h would seem to indicate that even with a representative level of noise added to the data, very good model recovery is possible. However, if the model parameterization is shifted laterally 25 km with respect to the synthetic model, the amount of recovery is significantly reduced. A shift-and-average inversion demonstrates that contrary to what would be inferred from examination of the resolution matrix or from conducting a synthetic spike inversion with velocity perturbations matching block boundaries (Figure 10h), discrete blocks of 100 x 100 x 50 km are not resolvable at depths of 100-400 km (Figure 10i). A synthetic test conducted using a single +5% spike of 150 x 150 km in the 200-260 layer showed that this was not an artifact caused by the choice of spike distribution. Synthetic anomalies spanning two layers (Figures 10i and 10j) are imaged. In an ideal case we would expect that the deduced wave speed in each block to be the volume-weighted average of the heterogeneity within it. However, this ideal is only achieved when illumination is uniform, as some rays must pass through any bit of heterogeneity in a block to signal its presence. In the case of non-uniform illumination, as is often the case in earthquake tomography, the result is difficult to predict. The above results suggest that an unfortunate choice of parameterization might mask heterogeneity on the scale length of a block, supporting the need for testing of multiple parameterizations. The recovery in the deepest layers of even the thin anomalies is surprising. It probably results from rays becoming more horizontal with depth and hence traveling a greater horizontal distance before crossing a layer boundary. However, it should be noted that these spike inversions might be giving an optimistic assessment of resolution in the deepest layers. At these depths the radius of the first Fresnel zone for a 1 Hz wave is ~75 km, but the effects of signal frequency are not included in the synthetic travel times generated by ray tracing. Certainly, there is little to be gained from using blocks much smaller than 100 km at these depths.
3.7. Interpretation and Evaluation of Velocity Model
Figures 11a-11m and Figures 12a-12d show horizontal and vertical slices through the block shift-and-average model. The ray segments close to the cross sections are also shown to provide an indication of the control of the travel time observations on the different portions of the model (please see our website at atlas.geo.cornell.edu (here) for depth slices of these ray segments). Regions that contain many crossing rays are better controlled than regions that contain a few subparallel rays.
The most robust feature imaged in the upper crustal layer (Layer 1) is the pronounced low-velocity
region seen in vicinity of the Strait of Gibraltar (Figure 11a)
both in Spain and Morocco. The low velocities are obtained despite the inclusion of significant station corrections
for stations in this region (see Figure 7). Inversions conducted
without station corrections yielded low-velocity perturbations of ~18% indicating that 8% is probably a minimum
value and that a significant portion of the anomaly may be being absorbed by the station corrections. These low-velocity
anomalies are significantly larger than the bootstrap standard deviation of ~3%. Low velocities in this region
were imaged by other tomography studies and refraction studies and, as mentioned earlier, may be caused by a combination
of low-velocity sediments and significant crustal deformation in the flysch regions. A-high velocity anomaly is
imaged in the eastern Rif to the SE of the Nekor fault (NF, Figure 2) in both crustal layers. This anomaly coincides with an aeromagnetic anomaly subparallel to the Nekor
fault that was interpreted by Michard et al. as indicating the presence of peridotite.
The most significant feature in the lower crustal layer (Layer 2) is the high-velocity anomaly
imaged off the Alboran coast of Spain in the vicinity of Malaga (Figure 11b). High velocities were also identified in this region by a refraction experiment . These authors reported
that ray trace modeling was unable to discriminate between a high-velocity peridotite body and a thin crust as
the source of this anomaly. The absence of stations in the Alboran Sea results in very poor control of its shallow
velocity structure with only limited rays passing through this layer beneath the Alboran Sea (See Figure
9). Gravity modeling and limited refraction profiling provide evidence of thinned (~20 km) crust,
but the underlying mantle velocities reported were highly variable ranging from 7.5 to 8.4 km/s on intersecting
profiles and therefore may be suspect. Low velocities are also imaged in the Rif and Betics.
Layer 3 from 30 to 45 km contains striking anomalies, consisting of a "horseshoe" of
anomalously low velocities (relative to an already low mantle velocity of 7.9 km/s) beneath the Betics, west Alboran
Basin (WAB) and Rif surrounding a region of high velocity beneath the Alboran Sea (Figure
11c). The low-velocity region may result from a mixture of two phenomena: (1) a crustal signature
resulting from crustal thicknesses of up to 38 km in places beneath the Betics and (2) the presence of anomalously
low-velocity mantle . The high-velocity region in the center of the Alboran Sea is surprising and difficult to
reconcile with our expectations of what to find beneath a region that underwent extension during the Miocene. One
possible geological explanation for this anomaly is that a significant thickness of thermal lithosphere has been
recovered since the termination of extension in the Middle Miocene. Although, Polyak et al. , based on heat
flow measurements, interpreted a lithospheric thickness in this region of only ~40 km. The central Alboran Sea
region contains few crustal earthquakes and effectively no stations. Another possible explanation is that it is
an artifact caused by assumption inherent in our approach. For example, we assume that a 1-D velocity model is
a good approximation for the first ray tracing and inversion iteration. The 1-D model places the Moho at depth
of 30 km. Rays connecting earthquakes on one coast of the Alboran with stations on the opposite coast, in reality,
are probably traveling at depths of ~15-20 km while they are traced through the model at depths of 30 km. This
extra segment of ray path before refracting along the Moho results in the theoretical arrival time being delayed
relative to the actual arrival times. However, as the rays are not passing through the 15 to 30-km layer above,
the velocities are increased in the 30 to 45-km layer in the least constrained region in the center of the Alboran
Sea. If this is true, then one must also consider the possibility that the surrounding lows might be enhanced by
the inversion attempting to compensate for the effect of these elevated velocities on longer ray paths. Another
possibility is that significant anisotropy may exist in the mantle and that inhomogeneous ray coverage might be
causing this anisotropy to be mapped into the isotropic velocity model. To check these hypotheses we conducted
1. Synthetics were generated using a mantle velocity of 7.5 km/s in the 15 to 30-km layer beneath
the Alboran. After adding noise, these data were inverted using the same starting 1-D velocity model. A significant
portion of the anomaly in the 15 to 30-km layer was mapped into the 30 to 45-km layer below, confirming that thinned
Alboran crust might be a viable explanation for this anomaly. However, although slight compensating lows were generated
beneath the Betics and Rif, they were minor compared to the very pronounced low-velocity anomalies imaged in the
2. As a further check on the independence of the Betic/Rif lows from this artifact, inversions
were conducted after Pn rays passing beneath the Alboran were removed. The African and Iberian phase data
were split into separate data sets and all local rays that intersected an E-W line passing through the Strait of
Gibraltar removed. These data sets were inverted separately using the same model parameterization. Although these
inversions used no ray paths traveling beneath the Alboran Sea, significant low-velocity anomalies were still imaged
beneath the Betics and Rif.
3. To investigate these anomalies further, as well as the possible influence of anisotropy, the
Pn arrivals extracted for various station-event separations were inverted for Pn velocity and anisotropy
using the code of Hearn and Ni . Inversions were conducted using varying amounts of Laplacian damping
on the isotropic and anisotropic velocity components to evaluate any trade-offs. Even when anisotropy was damped
very weakly, the low velocity beneath the Betics was imaged for all distance ranges. The Rif anomaly was also imaged
but found to be less robust than the Betic anomaly perhaps owing to the sparser ray coverage. Zero or positive
velocity perturbations were imaged beneath the Alboran Sea for all distance ranges and damping parameters. The
use of shorter ray paths did accentuate the high velocities beneath the Alboran Sea lending support to the thin
crust hypothesis; however, shorter paths might be less affected by a thin mantle lid, so this observation could
also be explained by the presence of a thin mantle lid beneath the Alboran Sea. Even with strong Laplacian smoothing
it was not possible to extend the imaged low velocities beneath the Alboran Sea, suggesting that uppermost mantle
velocities beneath the central Alboran may be higher than 7.9 km/s.
The 45 to 60-km layer (Layer 4) again shows low-velocity anomalies beneath the eastern Betics
and Rif (Figure 11d), but these lows are now separated by
a region of slightly elevated velocity beneath the western Betics and Rif. The low velocity beneath the Betics
conforms well to the surface outcrop of the Internal zone and its western boundary is marked by intermediate-depth
seismicity. Higher velocities are again imaged beneath the Alboran, but these anomalies are relatively small and
within one bootstrap error and so should not be considered robust.
The 60 to 100-km layer (Layer 5) is arguably the most important layer for evaluating any geodynamic
models, but also one of the most difficult to constrain as it marks the depth where ray coverage from local/regional
data is limited to only those from intermediate-depth earthquakes, and teleseismic rays are still limited to the
regions close to the stations (Figure 11e). The most robust
feature is the low-velocity zone along the southern coast of Spain that extends eastwards from the N-S line of
intermediate-depth seismicity. Despite the appearance that this anomaly extends beneath the Alboran Sea to the
African coast, this feature is predominantly controlled by teleseismic rays to the Spanish coastal stations that
do not sample the central Alboran. Beneath the western Betics and Rif a striking N-S high-velocity anomaly is imaged
to the west of the intermediate-depth seismicity. The bootstrap uncertainties indicate that the portion of this
anomaly beneath the flysch unit of the Betics is probably robust.
In layer 6 (100-150 km), a robust high-velocity anomaly is imaged beneath the WAB and west and
central Betics (Figure 11f). The SE extent of this anomaly
is not well defined owing to the absence of good ray coverage (Figure 9). Assuming that positive velocity anomalies may be taken to indicate a region of cooler mantle, rather
than a compositional difference, the most likely interpretation is that it indicates the presence of mantle lithosphere
in view of the coherent continuation of the anomaly to depth. The high-velocity anomaly beneath central Iberia
was not found to be stable using the bootstrap error test.
Layers 7, 8, and 9 (150-200-260-330 km) show a continuous progression of the high-velocity anomaly
to the southeast beneath the Alboran Sea (Figures 12g, 12h, and 12i). However, it should be noted that if the entire Alboran Sea is underlain by a high-velocity region at
these depths, this progression might in part be a reflection of the progressive improvement of ray coverage beneath
the Alboran Sea with depth.
The anomalies become relatively small and diffuse in layers 10 and 11 (330-410-490 km) (Figures 11j and 11k). However, significant concentrated anomalies are
imaged in layers 12 and 13 (490-570-650 km) (Figures 11l and 11m).
Any features imaged toward the base of the model must be interpreted with great care, as these regions are the
most likely to contain artifacts caused by deep mantle structure outside the model volume. One of the important
assumptions of an ACH-type approach is that any deviations from a spherically symmetric Earth outside the model
volume will be sampled equally by all the teleseismic rays from a single event. Masson and Trampert showed
that this assumption is sometimes not valid for a true heterogeneous Earth. Examination of a recent global tomography
model suggests that some heterogeneity may exist directly beneath the model volume but the size of imaged anomalies
is limited (variations of less than 1% from mean).
Cross sections AA', BB', CC', and DD' (Figure 12) provide a good summary of the major features seen in the velocity model and its relationship to seismicity. A low-velocity anomaly is imaged in the uppermost mantle along the SE coast of Spain. It is bounded on its western edge by intermediate depth seismicity. A high-velocity body appears to dip south and east from near lithospheric depths beneath Spain and the Strait of Gibraltar under the low-velocity region. This high-velocity anomaly is imaged to be continuous to depths of 350-400 km where it reduces in amplitude and becomes more diffuse. A pronounced high-velocity anomaly is again imaged at depths of 500-650 km near the base of the model in the vicinity of the deep earthquakes.
3.8. Synthetic Smearing Test and Comparison to Previous Velocity Models
The depth extent and volume of the high-velocity body in the upper mantle are important constraints for any geodynamic model for the region. Despite no explicit regularization between blocks, the smoothing inherent in both least squares and the shift-and-average methods will result in smearing of the imaged anomalies. To examine the effect of this smearing, synthetic travel times generated using a representative velocity model of the region (Figure 13) were inverted. Gaps were placed in key regions of the model (Layers 5, 10, and 11; Figures 13e, 13j, and 13k) to investigate whether these gaps could be recovered using the data distribution available.
Only limited smearing occurred into the 60 to 100-km layer (Figure
13e). A small positive anomaly was mapped into the region just to the west of the intermediate-depth
seismicity, in the vicinity of the imaged high velocity in the real model, suggesting that smearing may have played
some role in generating this anomaly. However, low velocities are not smeared into this layer from above, suggesting
that the low velocities beneath the southern coast of Spain in this layer are real. The existence of anomalous
low velocity mantle beneath the coastal Betics is further supported by the low Q values found for this region
by Canas et al. . Smearing clearly occurs into layers 10 and 11, indicating that if such a gap does exist
in reality at these depths, the inversion is probably biased toward filling it in.
Most previous studies have reported a relatively continuous high-velocity body from depths of
200 to 650 km . Examination of Blanco and Spakman's depth slices shows that similar to our model, a high velocity
is also imaged beneath the western Betics, Gibraltar, and Rif in their 70 to 120-km layer before moving east in
deeper layers to join with the deeper anomaly at 200 km. The studies by Plomerová et al. and Seber
et al. image a high-velocity body beneath the Alboran and southern Spain but owing to bottom layers beginning
at 200 and 150 km, respectively, they provide limited information on the depth extent of this body. Bijwaard
et al. show a section cut from their global model (BSE) in the Alboran region that also shows a pronounced
reduction in the size of the anomaly beneath 400 km (although this reduction in amplitude in the transition zone
may be a global feature of the BSE model (W. Spakman, personal communication, 1999). All these inversions use data
collected by permanent Spanish and Moroccan networks and so they possess the same potential for systematic smearing
errors. Establishing continuity of this body from 350 to 500 km is probably beyond the resolution of studies using
existing station distributions. However, existence of a single coherent body extending from < 200 km to 650
km is difficult to explain given the constraints on the relatively small amount of Tertiary convergence between
Africa and Iberia (~200 km). The hypothesis that the deep earthquakes occur in lithospheric bodies of limited size
is supported by the waveforms recorded in Spain. Buforn et al. suggest that puzzling arrivals from these
deep earthquakes may be caused by S to P conversions near the source and propose that they are occurring
in a body with a radius of only 26 km. A single body of this size would probably be unresolvable by our study.