3. Tomography Study

3.1. Data
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)

where

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 three tests:

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 actual model.

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.

Go Back