The upper-mantle transition zone beneath the Ibero-Maghrebian region as seen by teleseismic Pds phases

Abstract We investigate the upper-mantle transition zone thickness beneath the Ibero-Maghrebian region using teleseismic events recorded at 258 seismic stations from the IberArray of the Spanish TopoIberia project. For this purpose, the 410 km and 660 km depth discontinuities are mapped through the detection of P-to-s conversions. In order to add consistency and robustness to the detections, the final results are based on a joint analysis of receiver functions and two different cross-correlation functionals. Discontinuity depths are determined using time corrections obtained from a 3-D velocity model. This study presents the first high-resolution topography maps for the 410 and 660 discontinuities, which show lateral variations in the transition zone thickness beneath the study area. The results are discussed to add new constraints on temperature and composition to seismic velocity anomalies observed in the transition zone beneath the controversial Ibero-Maghrebian region. We obtained a thickened transition zone beneath the Alboran Sea and the north Balearic Sea suggesting low transition zone temperatures, and coinciding with the position of the high-velocity anomalies in the tomographic images related to the Betic-Alboran slab and the Alpine-Tethys remnant slab, respectively. A dominantly thinned transition zone beneath the Gulf of Cadiz and its surroundings suggest high-mantle temperatures. We also found evidences of the presence of a low-velocity layer atop the 410 discontinuity beneath the study area.

Tomography studies reveal a positive P-wave velocity anomaly beneath the Alboran Sea extending from the base of the crust across the entire upper mantle (e.g., Bezada and Humphreys, 2012;Faccenna et al., 2004;Monna et al., 2013;Piromallo and Morelli, 2003;Spakman, 1990;Spakman and Wortel, 2004;Villaseñor et al., 2003;Wortel and Spakman, 2000). This anomaly has an arcuate shape, but loses its geometry or resolution as the transition zone (TZ) is approached, at about 400 km depth. This heterogeneity has been explained through two different geodynamic models which also account for the formation of the Alboran Sea extensional basin. Different types of continental delamination or convective removal have been proposed (Calvert et al., 2000;Platt and Vissers, 1989;Seber et al., 1996) and, as opposite models, subduction of oceanic lithosphere with slab-tearing or partial slabdetachment (Bezada and Humphreys, 2012;Faccenna et al., 2004;Garcia-Castellanos and Villaseñor, 2011;Gutscher et al., 2002;Wortel and Spakman, 2000).
Independently from the tomography studies, SKS splitting (Buontempo et al., 2008;Díaz et al., 2010) and P-wave dispersion analyses (Bokelmann and Maufroy, 2007) are consistent with the presence of a subducted oceanic lithosphere (Bokelmann et al., 2011). Nevertheless, there is still room for the different interpretations. Additionally, recent tomographic images have revealed a negative P-wave velocity anomaly beneath the Gulf of Cadiz and south Portugal (Monna et al., 2013) extending to the upper-mantle transition zone, which might be related to high mantle temperatures and to the origin of the sublithospheric magma source responsible for the anorogenic magmatism in the Mediterranean.
Indirect evidences on upper-mantle temperature, composition, position, and vertical extension of the P-wave velocity anomalies would help their interpretation and can be achieved using seismic analysis. The study of the 410 km and 660 km depth discontinuities (or TZ discontinuities) is probably one of the best approaches, since these discontinuities are globally observed mineral phase transitions, which, as function of composition, respond with depth and thickness variations to temperature anomalies (e.g., Helffrich, 2000). Their study requires the detection and identification of seismic body waves which directly interact with the discontinuity through a reflection and/or wave type conversion (from P to S, or S to P). Here, we focus on P-to-s conversions, whose detection and extraction are based on their coherence, slowness, travel time, and polarity (e.g., Bonatto et al., 2013). In contrast to Bonatto et al. (2013), who studied the TZ at the Nubian-Eurasian plate boundary (Fig. 1a), we now investigate the transition zone beneath entire Iberia using a refined approach with about 6 times more data.
We present maps of the 410 km and 660 km depth discontinuities (hereafter referred to as 410 and 660 respectively) and the TZ thickness (TZT) beneath the Ibero-Maghrebian region, and we discuss the results in relation with the tomographic images from Villaseñor et al. (2003). Our study of the TZ discontinuities is based on the extraordinary dataset obtained from the TopoIberia seismic network (Díaz et al., 2009). A correlation between the anorogenic magmatism and the 410 and 660 topography is also investigated. Additionally, we study the presence of low-seismic-velocity zones atop the 410-km discontinuity, which are probably related to the presence of a melt layer caused by an increased water concentration in the TZ (Revenaugh and Sipkin, 1994;Schmerr and Garnero, 2007). Our analysis extends previous RF studies of the TZ discontinuities mainly beneath south Iberia (Bonatto et al., 2013;Dündar et al., 2011;van der Meijde et al., 2005), and the new data volume permits to resolve new TZ topography.

410-km and 660-km discontinuities
The 410 and 660 discontinuities define the limits of the TZ. In an upper mantle of pyrolitic composition (60% olivine with additional garnet and pyroxene), the 410 is the result of olivine to wadsleyite  Bird (2003). The yellow triangles mark the locations of the active anorogenic magmatism (alkaline magmatism of sub-lithospheric origin) (e.g., fig. 2 in Carminati et al. (2012) b) Location of seismic stations from the TopoIberia seismic network (red triangles) and from other collaborating institutions (yellow triangles). Fig. 2. Distribution of earthquakes used in this study. White lines correspond to circles located every 30 ∘ and centred at the yellow triangle, which denotes the location of the TopoIberia seismic network (IberArray). transition, while the 660 is the dissociation of ringwoodite into perovskite + magnesiowüstite. Temperature anomalies in the mantle move the phase changes to different pressures (depths) according to the Clapeyron slopes (Bina and Helffrich, 1994). Due to the opposite sign in the Clapeyron slopes of the phase changes responsible for each discontinuity, the 410 and 660 depth changes are anti-correlated as response to a thermal anomaly. While the 410 becomes shallower in colder regions and deeper in hotter ones, the 660 behaviour is opposite. As a consequence, the TZ becomes thicker near subducted slabs and thinner beneath plumes or high-temperature regions due to small-scale mantle convection (Collier et al., 2001;Helffrich, 2000;Lawrence and Shearer, 2006;Vidale and Benz, 1992). However, the presence of other transforming or non-transforming mantle components can change the characteristics of the phase transitions and makes the interpretation of the corresponding seismic discontinuities more difficult (Cao et al., 2011;Deuss et al., 2006;Stixrude, 1997;Tauzin et al., 2008;Thomas and Billen, 2009;Wang et al., 2004;Weidner and Wang, 1998). In particular, in normal to hot upper mantle and in the presence of Al within garnet, the garnet-to-perovskite transition is expected near the bottom of the TZ with a positive Clapeyron slope (Akaogi et al., 2002;Wang et al., 2004;Weidner and Wang, 1998). As a consequence, the garnetrelated 660 discontinuity deflects downward in a hotter mantle, contrary to the olivine-related 660. In addition, the depth of the garnet-related phase transition is strongly dependent on the Al content of garnet, causing variations in the discontinuity depth by more than 50 km (Weidner and Wang, 1998).

Data and method
The data used are from 258 stations shown in Fig. 1b, which belong to the TopoIberia (Díaz et al., 2009) seismic network (red triangles), to the WILAS project, and other collaborating institutions (yellow triangles). The data set is composed of 509 events from epicentral distances between 65°and 95°and a magnitude mb between 5.5 and 7.5 for a time period of 5 years (from January 2007 to December 2011). Fig. 2 shows the worldwide distribution of these earthquakes. In this study, we use events with epicentral distances between 65°and 95°, since the majority of events are at these distances. Besides, the distance range used avoids interference with phases such as PP and PcP. As the stations were not simultaneously active, not all the events were recorded in all the stations. The selection of events is based on the SNR in the vertical (Z) component and on the subsequent visual check of the radial (R) component.
The 410 and 660 are commonly studied through the detection of P-to-s conversions (Pds, where d is the depth of the discontinuity) ( Fig. 3a) on teleseismic events (e.g., Eagar et al., 2010;Li and Yuan, 2003;Paulssen, 1985;Vinnik, 1977). Differential travel times between P-to-s conversions and P provide a measure of the one-way S travel time between the surface and the discontinuity. Therefore, using an adequate velocity model, this one-way travel time can be translated to discontinuity depth. P-to-s conversions from mantle discontinuities arrive in the P-wave coda of the R component of teleseismic events. These conversions arrive together with multiply reflected and scattered waves, and are difficult to observe directly on seismograms mainly because of their low amplitude and abundance of other signals. In addition, P-to-s conversions are expected to be coherent with the waveform of the first arrival for conversion at discontinuities which are thinner than one half of the wavelength (Bostock, 1999;Paulssen, 1988;Richards, 1972).
In order to extract the P-to-s conversions by the waveform similarity, three different approaches are used following Bonatto et al. (2013). We use deconvolution of the Z component from the R in the frequency domain, which is known as the receiver function (RF) technique (Ammon, 1991;Langston, 1979;Phinney, 1964;Vinnik, 1977). Deconvolution of the Z component from the R is equivalent to spectral division in the frequency domain. To stabilize the spectral division, we use the waterlevel deconvolution (Clayton and Wiggins, 1976). Additionally, the waveform similarity between the P-to-s conversions and the P phase is inspected with two independent cross-correlation techniques: the classical cross-correlation geometrically normalized (CCGN) and the phase cross-correlation (PCC) presented by Schimmel (1999). These techniques are based on different strategies. In analogy to the classical crosscorrelation, the PCC measures the waveform similarity between two signals as function of lag time. The PCC is based on the instantaneous phase similarity of the corresponding analytic traces, is amplitude unbiased, and is more sensitive to waveform coherence than CCGN (Schimmel, 1999;Schimmel et al., 2011) and, therefore, well suited for the detection of coherent weak-amplitude signals. CCGN is based on the sum of signal Fig. 3. a) Example ray-path for a P-to-s conversion at the TZ discontinuities at 410 km or 660 km depth. P wave is shown as continuous line, S wave as discontinuous line. b) Phase weighted stack of 165 radial cross-correlograms (PCC and CCGN) and receiver functions (RFs) in the relative time-slowness domain, for a reference distance of 80 ∘ . amplitude products and is therefore less sensitive to waveform coherence. The decreased sensitivity may favor signal detection when there is less waveform similarity due to waveform distortion. Both crosscorrelations are performed between the R component and a portion of the Z component which contains the P phase and part of its coda. Thus, the Z component of each trace is windowed from the onset of the P phase to 100 s after the first arrival. The same window is used to compute the RFs. The R and Z components are band-pass-filtered between 0.03 Hz and 0.2 Hz prior deconvolution (or cross-correlation).
Independent quality controls are applied over individual RFs and PCC-and CCGN-correlograms. These controls are based on the SNR, which is defined as the rms ratio of the absolute amplitudes in a 30 s window before and after the P phase arrival in the R component of the corresponding PCC, CCGN, or RF. After the quality control, we obtain 8819 RFs, 8085 R-correlograms with CCGN, and 8027 with PCC. Finally, the P-to-s conversions are detected on stacked RFs and cross-correlograms (PCC and CCGN). We use the phase weighted stack (PWS) (Schimmel and Paulssen, 1997) in the relative time-slowness domain, where the reference is the P phase, as shown in Fig. 3b for a subsidiary data set.

Stacking and binning of data
In order to stack, the data is grouped into bins of common piercing point (CPP) areas. Fig. 4a shows the piercing points of P-to-s conversions at a discontinuity at 510 km depth obtained with the 8085 event-station pairs used to compute the CCGN-correlograms. This figure illustrates the approximate region in the upper mantle sampled by P410s and P660s phases. For a reference distance of 80°, the conversion point distance for P510s-P660s is 40 km and for P510s-P410s is 30 km. The lateral resolutions of P410s and P660s as defined by their corresponding Fresnel zones at 5 s period (0.2 Hz) are circle-like areas with radii of 120 km and 180 km, which is more than twice the distance to the piercing points mapped at 510 km depth. Here the Fresnel zone is defined as the area where the travel time differs by less than a quarter period. Consequently, the Fresnel zone at 0.2 Hz is a lower bound since 0.2 Hz is the highest frequency used in our data.
The piercing point density depends on the station and event distribution and is variable over the entire study region. This implies a variable discontinuity visibility, which we account for using adaptive bin sizes depending on the local number of piercing points. The bin sizes vary between 0.8°, 1°, and 1.6°width in latitude and longitude and are located every 0.25°in both directions (see example of the bin size and overlap in Fig. 4a). The bin sizes are on the order of the Fresnel zone of P410s and P660s at a period of 5 s. Fig. 4b shows the grid for the CCGN-correlograms, where the different colours discriminate between bins of different sizes. The grid is composed of 2749 bins for PCC, 2751 for CCGN, and 2766 for RF. The number of piercing points per bin which equals the number of CCGNcorrelograms to be stacked for each bin is shown in Fig. 4c. This figure shows the available data density and areas where we expect a higher robustness of results due to the higher discretization into independent data bins. The best data coverage is in the southern Iberian region and in the Pyrenees. Small differences in the data coverage have been observed with the different techniques (CCGN, PCC, and RF), which, however, do not have a meaningful influence on the results.

Detection of Pds phases and quality criteria
For an unambiguous detection of P410s and P660s, the data are stacked per bin using the PWS in the relative time-slowness domain using a reference distance of 80°(as shown in Fig. 3b). The reference distance of 80°for the slant stacks is chosen to account for the distribution of epicentral distances which range from 65°to 95°. Bins with less than 12 traces are not considered for this analysis. The P410s and P660s phases are extracted by means of their coherence (or amplitude for RFs) and their expected slowness, travel time, and polarity. The signals are identified as positive maxima with amplitude larger than twice the absolute mean amplitude of the stack in the time interval 30-80 s. Stacks without signals inside the time intervals 42.8 s ± 5 s and 65.8 s ± 5 s and the relative slowness interval of ± 0.4 s/°have been discarded. These time intervals are centred at the P410s and P660s relative times in AK135 (Kennett et al., 1995) for a reference distance of 80°.
A bootstrapping technique (Efron and Tibshitani, 1986) of 21 repetitions is used to estimate the uncertainties in our relative travel time measurements, tP410s-tP and tP660s-tP. This is done for each bin and processing technique (PCC, CCGN, and RF) to assess the robustness of the individual detections per bin. Signals with bootstrap uncertainty larger than 1.5 s (~15 km) are discarded. Most of the stacks that pass the quality controls show only one clear signal near P410s and P660s (e.g., Fig. S1c in the supplementary material). However, some stacks present more than one signal near the reference time of P410s and P660s (e.g., CCGNs and RFs stacks in Fig. S1a or RFs stack in Fig. S1b). Therefore, we perform a visual inspection of the time-slowness stacks Fig. 4. a) P510s piercing points (grey crosses) in the study area. The squares show the bin overlap and size. b) Bin size for CCGN-correlograms: blue = 0.8 ∘ width; red = 1 ∘ width; yellow = 1.6 ∘ width. c) Data density for the bins shown in (b). The best data coverage is obtained in south Iberia and in the Pyrenees. of the three techniques simultaneously in order to decide which of the detected signals is actually the P410s and/or P660s phases, or to discard data when a clear identification of phases is not possible. Note that after the quality controls, we may remain with bins whose measures are based on only one or two techniques. In Fig. S2 (supplementary material), we show additional examples of signal detection which were selected based on the clear visual detection of both phases.

Integrated detections from the different techniques
Whenever possible, we merge the independent P410s and P660s detections from each bootstrap repetition for each of the three techniques. This way, new mean and standard deviation values are obtained for the relative travel times tP410s-tP and tP660s-tP at each latitudelongitude bin. Only travel times with bootstrap uncertainty smaller than 1 s (~10 km) are considered for the final analysis. Finally, 923 and 1627, respectively are the number of bins with robust P410s-and P660s-detections.

Time corrections and depth conversions
We use the P-wave tomography model by Villaseñor et al. (2003) (an updated tomography model from Wortel and Spakman (2000) due to added events for Spain) to account for volumetric seismic velocity anomalies beneath the study area and to correct the estimated relative travel times before depth conversion as in Bonatto et al. (2013), Dueker and Sheehan (1997), and van der Meijde et al. (2005). We compute average 1-D P-wave velocity profiles for 2°× 2°bins centred at each CPP bin. In analogy, S-wave velocity profiles are derived from the P-velocity anomalies, δv p , by employing a constant factor δv s / δv p = 1.5. This factor typically ranges between 1.5 and 2 in the upper mantle (Ritsema and Van Heijst, 2002), depending on the type of anomaly (thermal and/or compositional). The corresponding time corrections are computed with respect to AK135 and range between −1.2 to 0.3 s for the P410s phases and − 1.2 to 0.6 s for the P660s phases. Changing the applied constant velocity perturbation ratio δv s / δv p = 1.5 to 2 leads to high end time corrections in the order of − 1.8 s to 0.9 s, which translates to depth variations of about −6 km to 3 kmwith respect to using δv s / δv p = 1.5. Finally, the corrected time values are converted to depth using AK135 and an 80°reference epicentral distance. A more accurate depth correction is not performed due to the highly variable resolution of the uppermantle velocity structure beneath Iberia and Morocco.
The corrected relative travel time values for P410s and P660s at each CPP bin and the corresponding standard errors are shown in Fig. 5. The corrected values have been plotted relative to the time values in AK135; that is, the plotted values are (tPds − tP) Estimated − (tPds − tP) AK135 (d = 410,660). For P660s, blue areas correspond to converted phases arriving after the reference value and red areas to converted phases arriving before the reference. In order to enhance the opposite effect of temperature on the topography of discontinuities a reversed colour bar is used for P410s. Note that there are several locations without P410s time values (in grey). The origin of these gaps is either due to the lack of P410s-detections, or due to detections of incoherent lowamplitude signals which do not satisfy the employed quality criteria. This lack of detection might be due to structural complexities that cause defocusing or loss of coherence and destructive interference with other coda waves.

Results
The distribution of our estimated depth values for the 410 and 660 discontinuities before time corrections is shown in Fig. 6a and b with grey-filled bars. The corresponding average values of the apparent discontinuity depths are 416 km and 669 km (black line). After time correction, the average depth values move to 412 km and 663 km (red lines), which are closer to the reference depths of 410 km and 660 km (grey lines). The upper labels in Fig. 6a and b show the predicted temperature variations using Clapeyron slopes of + 3.1 MPaK − 1 for the olivine-wadsleyite transition (410) and − 2.6 MPaK − 1 for the post-spinel transition (660) (Akaogi et al., 2007).

410 and 660 maps
The interpolated discontinuity depths are presented in Fig. 7. We use the nearest-neighbour interpolation algorithm from the GMT plotting tool (Wessel and Smith, 1991) to assign an average value to each node. The average value is computed as a weighted mean of the nearest point from each sector inside the search disk with a radius of 150 km. This final smoothing is like averaging to identify the main and more robust topography trends. The yellow triangles in the 410 map in Fig. 7 show the locations of the active anorogenic magmatism (e.g., fig. 2 in Carminati et al., 2012), of sub-lithospheric origin and occurring after the volcanic episodes associated with the subduction (orogenic magmatism) (Lustrino and Wilson, 2007;Lustrino et al., 2011). This activity is mainly Pliocene to recent in age. In this context, active means that these volcanoes erupted recently (from a geological point of view) and could still erupt. Note that they are all located at the border of areas with a depressed 410. The yellow star in the 660 map in Fig. 7 shows the epicentre of a nest of deep earthquakes beneath Granada (e.g., Buforn et al., 2004), which coincides with a transition from depressed to uplifted 660.

Transition zone thickness
Differential travel times tP660s − tP410s are mainly sensitive to 3-D heterogeneities within the mantle TZ because both phases, P410s and P660s, are nearly identically affected by lateral variations in the Earth's structure above 410 km. Besides, velocity anomalies within the TZ are expected to be smaller than in the uppermost part of the mantle. Therefore, the analysis of the TZT yields a more accurate result than the individual 410 and 660 depths. The TZT values are computed from the subtraction of the 410 and 660 depths (TZT = H660 − H410) after interpolation. Fig. 6c shows the distribution of apparent TZT (grey-filled bars) and the TZT after applying the travel time corrections to the P410s and P660s phases (red-unfilled bars). The averaged value of the apparent TZT is 252 km, slightly thicker than the AK135 value of 250 km. After time corrections the averaged TZT is 250 km. Both values are in good agreement with the global average value of Chevrot et al. (1999).   (Buforn et al., 2004). In order to enhance the opposite effect of temperature on the topography of discontinuities a reversed colour bar is used for P410s. Black lines a, b and c denote the location of the three profiles shown in Fig. 10. Black dots show the locations with detection of P410s and P660s.
The TZT is shown in Fig. 8. The green squares show the location of selected bins with robust 410 and 660 depth values, where the P410s and P660s phases have been detected with at least two techniques and with small uncertainty values. These bins are located in areas of special interest, where the TZ shows thickness values larger or smaller than the reference thickness in AK135. The estimated 410 and 660 depths and the TZT for the green squares in Fig. 8 are listed in Table 1.

Average TZT
In the central part of Iberia, the TZT shows variations around the mean of about ±10 km (A in Fig. 8). Also the 410 and 660 topographies show depth variations in central Spain smaller than ± 10 km (A in Fig. 7 a  and b).

Thinner TZ
The southwest coast of Iberia (C2 and D1 in Fig. 8 and Table 1) exhibits a thinner TZ of about 220 km. With the exception of the Spanish side of the Strait of Gibraltar, the thinner TZ in this region is not due to anti-correlated 410 and 660 depths. The 410 is, on average, deeper in the whole area, but the largest depressions are towards the east, beneath the Gulf of Cadiz and the Strait of Gibraltar. On the other hand, the 660 is, on average, deeper to the east and shallower to the west of the region. In this area, the tomographic images from Monna et al. (2013) show a negative P-wave velocity anomaly extending to the upper-mantle transition zone (from sub-lithospheric depths down to the bottom of the model at about 660 km depth), interpreted as hot mantle.

Thicker TZ
The southern and eastern coasts of Spain are characterized by a thicker TZ with maximum thickness of about 300 km beneath the Alboran Sea (B in Fig. 8 and Table 1) and about 280 km beneath the Pyrenees and north Balearic Sea area (F1 in Fig. 8 and Table 1). In these areas, the tomographic images (Spakman and Wortel, 2004;Villaseñor et al., 2003;Wortel and Spakman, 2000) show positive P-wave velocity anomalies related to the Betic-Alboran slab and to the Alpine-Tethys remnant slab, which seems to be stagnant at the base of the TZ (e.g., feature A in profiles k to l in fig. 2 A2.2 in Spakman and Wortel (2004)). Beneath the Alboran Sea, the TZ is thick due to an elevated 410 and a depressed 660. The broad extension of the thicker TZ beneath the Alboran Sea is controlled by the extended 660 depression, while the 410 is locally elevated in the areas where the TZ reaches the largest thickness value. On the other hand, beneath the eastern Pyrenees and the northern Balearic Sea both discontinuities are depressed. A thick TZ is also visible beneath the northwest of Morocco (D2 in Fig. 8 and Table 1), primarily due to an up to 20 km depressed 660 discontinuity.
Another interesting result is that the deep earthquakes beneath Granada coincide with a transition from thinner to thicker TZ. In this area, our observations present a large P410s-detection gap and a smaller P660s gap. Therefore, the TZ is obtained exclusively from interpolated 410 and 660 values.
In the north and northwest coast of Iberia the 410 and 660 depths are correlated. Both discontinuities are deeper beneath the Bay of Biscay (E in Fig. 7) and shallower beneath the northwest coast of Iberia (feature C1 in Fig. 7). In both areas, the TZT is close to the reference value of 250 km. It is probable that the 410 and 660 depth correlation is due to underestimated time corrections on P410s and P660s travel times. This is expected due to the lower resolution of the tomographic images in these areas. For this reason, these areas are not further investigated.

Negative amplitude precursory signals
We also observed negative amplitude signals before P410s and P660s. Since the RFs hold the relative amplitude information, we investigate these additional features in the cross-sections of RFs. Representative examples of these negative amplitude precursory signals are shown in the four cross-section of stacked RFs in Fig. S3 (in the supplementary  material).
The negative amplitude signal which accompanies the P410s phases is denoted as Pws (marked with red lines in Fig. S3). In some areas, the amplitudes of these particular signals are comparable to that of the P410s phases (e.g., a in Fig. S3). In other places, the only signal which is clearly seen is the negative one, while the P410s amplitudes are significantly diminished, e.g, beneath the western Alboran Sea (c in Fig. S3). Akin precursors have also been observed by other researchers (e.g., Eagar Table 1. Black dots show locations with simultaneous detection of P410s and P660s. The TZT is obtained from the subtraction of the 410 and 660 depths after interpolation using the whole data set. in other regions, and have been interpreted as caused by a melt layer (of low velocity) which is often explained through an increased water content in the TZ (Bercovici and Karato, 2003;Tauzin et al., 2010). Therefore, we attribute the negative precursory signal Pws to a conversion from a first order decrease in velocity atop the 410 discontinuity. By waveform modelling, Pino and Helmberger (1997) have developed a one-dimensional model for the upper mantle compressional velocity structure in the West Mediterranean Basin. In consistency with our preferred explanation, they observed pronounced low velocity gradients between 200 km and 368 km depth suggesting the presence of melt or fluids associated with subduction. The negative amplitude arrival before the P660s phases is denoted as X1 in Fig. S3. Up to our knowledge, there exists no known mineralogical phase change that would produce the negative signal before P660s. Other studies have also observed similar negative impedance features and have interpreted them as accumulated oceanic crust at the base of the TZ (Eagar et al., 2010;Shen and Blum, 2003;Shen et al., 2008Shen et al., , 2014. This signal is not used in the following analysis, and therefore, we do not further explore its origin.

Relation with previous works
The large data coverage of IberArray (Díaz et al., 2009) has permitted us to study in detail and with high resolution a previously under-sampled portion of the TZ discontinuities beneath the Ibero-Maghrebian region. Our results locate the 410 and 660 discontinuities within the expected depth range as obtained by global studies (Chevrot et al., 1999;Lawrence and Shearer, 2006;Shearer, 1991;Tauzin et al., 2008;Vidale and Benz, 1992), and are further in good agreement with a large scale P-RF study that covered the entire Mediterranean region (van der Meijde et al., 2005). In van der Meijde et al. (2005) two stations are located in central Spain and north Morocco, showing the same TZT variations as observed in our study: small thickness variations beneath central Spain and a thicker TZ beneath Morocco. The P410s-detection gap and the TZ thickening confirms the observations by Bonatto et al. (2013), which have been based on less data and a larger bin size and spacing. The observation of the P410s-detection gap confirms also the study of Paulssen (1988), van der Meijde et al. (2003), Morais et al. (2015), which also show intermittency for such an interface beneath Iberia. The thickening of the TZ towards Africa seems to be a robust feature; we have estimated a TZT variation from 214 km to 270 km in the Strait of Gibraltar (D1 and D2 in Fig. 8).

Integrated detections
The joint usage of the three independent approaches permits us to stabilize the detections against measurement errors, to use the measurement variability as a robustness indicator, and to bridge observation gaps (Bonatto et al., 2013). The Venn diagram shown in Fig. 9 illustrates the distribution of the number of bins with robust detections of P410s (left panel) and P660s (right panel) with the different techniques, PCC, CCGN, and RF. The diagram shows that we have obtained 369 bins with P410s-detections for PCC, 564 bins for CCGN, and 689 bins for RF. The number of bins with P660s-detections is larger, we obtain 930 detections with PCC, 1235 with CCGN, and 1284 with RF. When the results from the different techniques are considered, 923 and 1627, respectively are the number of bins with robust P410s-and P660s-detections , which are larger than the number of detections with the individual techniques. Note that there are always a certain number of bins which show robust detections with only one technique (numbers outside intersected areas). This fact enhances the importances of using different approaches to increase the number of travel times used to compute the 410 and 660 topography maps.

Quantitative analysis of time corrections and discontinuity depths
In a mantle dominated by the olivine phase transitions, the depths of the 410 and 660 are anti-correlated in both cold and warm regions. Nevertheless, the mantle structure leads to large positive correlation of the observed depth values (and relative travel times), e.g., +0.94 in Chevrot et al. (1999), +0.84 in Gao and Liu (2014), and +0.6 in Tauzin and Ricard (2014). Ideally, a negative or near zero correlation value is expected after correction, but in practice it is observed that velocity corrections reduce the correlation to smaller positive values (e.g., Gao and Liu, 2014;Tauzin, 2009). This residual positive correlation is partly due to the limitations of the velocity models to reproduce the real structure of the Earth, due to the assumption that depth changes are controlled only by temperature changes, and to a lesser extent due to measurement errors. Fig. S4 shows the relation between H660 (depth of the 660) and H410 (depth of the 410) before (a) and after (b) applying time corrections. These figures were obtained with the results from the CPP bins with simultaneous detection of both phases, P410s and P660s. The grey line represents perfectly correlated topographies. After time correction, the correlation between the corrected depth values reduces from 0.34 to 0.23 and the mean value of (H410, H660) denoted with the red point moves toward the theoretically expected value (grey square). The decrease of correlation implies that the influence of mantle velocity structure is minimized after correction. A small correlation value before correction, as the one we obtained, may be related with a heterogeneous distribution of (H410s, H660) values. This is expected in studies covering smaller areas. For example, Tauzin and Ricard (2014) studied an area beneath west USA which is contained in the study area of Gao and Liu (2014) and found a smaller correlation value than in Gao and Liu (2014).
It is widely assumed that velocity perturbations in the tomography models represent mainly temperature changes. Therefore, in an olivinedominated phase-transition system, a negative correlation is expected between the corrected 410 depth variations (dH410 = H410 − 410km) and the velocity anomalies and a positive correlation are expected for the 660 depth variations (dH660 = H660 − 660km). Before correction, the velocity structure tends to amplify the 410 topography and to attenuate or reverse the topography of the 660. As a result, the correlation value before correction is expected to be large and negative for the 410 and near zero or negative for the 660 (e.g., Tauzin and Ricard, 2014, fig. 1). Fig. S5a shows the estimated dH410 after correction versus the average P-wave velocity anomaly above the TZ from Villaseñor et al. (2003). The correlation coefficient is R = 0.05, while before correction it was R = −0.2. This correlation decrease means that the influence of the velocity anomalies above the TZ on the 410 depth has been reduced through the time correction. Similar conclusions are derived from Fig. S5b, which shows the estimated dH660 after correction versus the average P-wave velocity anomaly in the top 700 km. Nevertheless, the values in Fig. S5a and b before correction are very low compared to other studies (Chevrot et al., 1999;Gao and Liu, 2014;Tauzin and Ricard, 2014). These poor correlation values probably indicate that depth changes are not controlled primarily by temperature changes and/or that the velocity structure in the upper-mantle and the TZ are anti-correlated with anomalies of the same order but opposite sign.
Finally, in order to evaluate the influence of the computed corrections on the final relative time values, we also show the relation between the time corrections or predicted residuals (tPds AK135 − tPds tomography ) and estimated residuals (tPds AK135 − tPds obs,corrected ) for P410s (Fig. S6a) and for P660s (Fig. S6b). The absence of correlation between predicted and observed time residuals strengthens that we are not interpreting topography produced by the time corrections.

Alboran Sea and surrounding areas: thicker TZ
The elevated 410 and depressed 660 in the Alboran Sea lead to a 300 km thick TZ, which is consistent with the position of the positive P-wave velocity anomaly of the Betic-Alboran slab (e.g. Bezada et al., 2013;Villaseñor et al., 2003;Wortel and Spakman, 2000). Under the assumptions of a 250 km TZ thickness for reference, a pyrolitic mantle composition, and Clapeyron slopes of + 3.1 MPaK − 1 for the olivinewadsleyite transition (410) and −2.6 MPaK −1 for the post-spinel transition (660) (Akaogi et al., 2007), the TZ thickening of 50 km beneath the Alboran Sea can be translated into an approximate temperature decrease of 292 K (Helffrich, 2000). A 10 % change of the Clapeyron slopes leads to a variation of about 10 K. The Betic-Alboran slab is probably of oceanic origin (hydrated oceanic lithosphere) (Bokelmann et al., 2011), then, the presence of signal Pws in this area (b and c in Fig. S3 between about latitudes 35°N and 38°N) may be related to an increased water content in the TZ. Hence, part of the 410 uplift and 660 depression could be due to a water content increase in wadsleyite and ringwoodite (Litasov et al., 2005). Therefore, the temperature anomaly of − 292 K should be considered as a higher bound.
The composed tomography-discontinuity profile at 3°W longitude ( Fig. 10a) clearly shows that the uplifted 410 and depressed 660 beneath south Spain (B in Fig. 7) are vertically anti-correlated, as illustrated in model 1 in Fig. 11. This implies that the thermal anomaly is likely to extend vertically in depth through the entire TZ, which is consistent with the steep Betic-Alboran slab (Bezada et al., 2013).

The area of deep earthquakes beneath Granada
The east-west profile in Fig. 10b, centred at latitude 37°N, shows that the depression of the 660 includes an area of isolated deep earthquakes beneath Granada (37°N 3.7°W approximately) with similar focal mechanism, magnitudes (M W ) from 4 to 7, and a recurrence time of about 20 years (Bezada and Humphreys, 2012;Buforn et al., 2004Buforn et al., , 2011. Furthermore, Fig. 8 shows that the location of the nest of deep earthquakes (yellow star) corresponds to a TZT change from thin to thick, which may be interpreted as a temperature change. The thinner TZ in this area is due to anti-correlated 410 and 660 depth changes. Assuming olivine-related TZ discontinuities and using the values in Table 1, we estimate a net temperature variation of about 502 K (from about − 292 K inside the cold slab to about + 210 K at the western end of the slab) over a distance of less than 250 km inside the TZ. That is, these earthquakes occur in an area of high temperature gradient.

Pyrenees and northeast coast of Spain
The depressed 660 in this area (F1 in Fig. 7) coincides with the location of the positive P-wave anomaly related with the Alpine-Tethys remnant slab (e.g., feature A in fig. 2.4 in Spakman and Wortel (2004)). Thus, we interpret the 20-30 km downward deflection of the 660 as caused by the cold slab, which shifts the post-spinel transition to greater depths. We infer a negative thermal anomaly of about 256-384 K at the 660, which we consider as a higher bound consistent with the assumption that part of the 660 depression is caused by the presence of water (Litasov et al., 2005), perhaps from dehydration of the slab. In this area, the 410 is also downward deflected by about 20 km (F1 in Fig. 7) in consistence with the presence of low P-velocity above the TZ in this area (see Fig. 10c). Therefore, the 410 downward deflection may be related to high temperature above the stagnant slab. Nevertheless, as will be discussed in Section 6.9, we cannot discard a compositional explanation for the origin of the 410 downward deflection. Additionally, strong and clear Pws phases are observed in this area (a in Fig. S3), suggesting water-induced melting atop the 410 (e.g., Tauzin et al., 2010).
The Alpine-Tethys remnant slab is of oceanic origin and probably contains a high proportion of water inside. It is probable that smallscale convections, triggered by slab dehydration (Richard and Bercovici, 2009), occur above this stagnant slab and provide the water to hydrate the olivine above the TZ. In this context, the presence of hydrated wadsleyite may explain the reversed-polarity signal before P410s. The preferred model for this region is model 2 in Fig. 11, where the grey arrows mark the locations of the small-scale convection which transports the water responsible for the water-induced melting atop the 410.

Thicker TZ beneath the western part of the Moroccan region
The maximum observed TZT beneath Morocco (275 km) is 25 km thicker than global averages of 250 km (Chevrot et al., 1999;Flanagan and Shearer, 1998;Lawrence and Shearer, 2006). This thickness increase is mostly due to a deeper 660, which shifts downward by as much as 20 km (D2 in Fig. 7 and in Table 1), while the 410 shows less P410s detections (Fig. 5) and stays on average at its expected nominal depth of 410 km. Under the assumption of an olivine-dominating composition, the TZ thickening of 25 km beneath Morocco can be translated into an approximate temperature decrease of 146 K. However, there is no corresponding positive velocity anomaly in the tomographic image by Spakman and Wortel (2004), which could strengthen this hypothesis.
A recent tomography study (Bezada et al., 2013), which includes more stations in Morocco, has revealed a high-velocity body at a depth of about 400-500 km beneath the Middle Atlas, which is probably a portion of delaminated lithosphere. This high velocity anomaly would only decrease the travel time of the P660s phase, which, without the correct travel time corrections, would lead to apparent discontinuity uplift. However, we observe the opposite, a deeper 660. Thus, a plausible cause for the thickened TZ can be an Al-rich garnet-related 660 (Wang et al., 2004;Weidner and Wang, 1998) as proposed in Bonatto et al. (2013). In the absence of a hot thermal anomaly, as expected from the tomographic images, the depth variations of the garnet-related 660 beneath Morocco can likely be attributed to variable proportions of Al content within garnet (e.g., Thomas and Billen, 2009).

Thinner TZ beneath Strait of Gibraltar, Gulf of Cadiz and south of Portugal
In this area, the P-wave anomalies are slightly negative (Spakman and Wortel, 2004;Villaseñor et al., 2003) indicating velocities lower than the reference model, which suggests that the mantle is warmer than average. Recently, Monna et al. (2013) published a tomography study which also covers the Atlantic domain in the southwest Iberian margin. This study has a higher resolution than previous studies because it includes data from OBS stations in the Gulf of Cadiz and the Atlantic Sea. Below the Atlantic domain, they image a wide area of negative P-wave velocity anomalies down to the base of the TZ, which is attributed to the presence of a hot upper mantle and which is perhaps related with the late Mesozoic early Cenozoic alkaline igneous activity of sub-lithospheric origin (e.g., Lustrino and Wilson, 2007;Grange et al., 2010, and references therein).
With the exception of the Spanish side of the Strait of Gibraltar, the thinner TZ in this region is not due to anti-correlated 410 and 660 depths. The 410 is, on average, deeper in the whole area, but the largest depressions are towards the east, beneath the Gulf of Cadiz and the Strait of Gibraltar. On the other hand, the 660 is, on average, deeper to the east and shallower to the west of the region, but it has small length-scale depth changes (note the depth variability of the 660 in the western side of Fig. 10b). The depth variability of the 660 may suggest the presence of multiple discontinuities, but, due to the used quality criteria, we have detected each time the best quality signal.
Our preferred model for this region is model 3 in Fig. 11, where a warm thermal anomaly crosses the TZ causing the depression of the 410. At the base of the TZ, the presence of multiple discontinuities cannot be ruled out.

410 depression and anorogenic magmatism
The recent anorogenic magmatism (b5 Myr) (Carminati et al., 2012), depicted with yellow triangles in Fig. 7, are all located at the border of areas with a depressed 410, while no spatial correlation is seen with the 660 topography. All these volcanic centres should possibly be viewed within the wider context of the Na-rich intra-plate magmatism (since late Cretaceous to present) characterizing the whole Euro-Mediterranean region and extending as well to the eastern Atlantic realm (e.g., Lustrino and Wilson, 2007;Downes, 1992, 2006). Interestingly, the distribution of all alkaline volcanic centres which have been active in the last 30 Myr has a first order correlation with low P-velocity anomaly regions averaged over the 100-400 km depth range in the upper mantle (Piromallo et al., 2008, fig. 7a). Furthermore, the bulk of the Euro-Mediterranean alkaline volcanic fields active in the last 30 Myr is located at the northern edge of the widespread high velocity anomaly in the TZ (Piromallo et al., 2008, Fig. 7b), where cold material is accumulated since the Eocene by different subductions. Piromallo et al. (2008) proposed that the upward return-flow driven by subductions and related extensional processes likely caused intermittent alkaline activity. In this context, the signature of previously contaminated mantle (polluted at the time of a larger-scale upwelling from the lower mantle) will emerge at different times and places, particularly where local tectonics favor its rise (e.g., rift areas, back-arc regions, slab windows). Lithospheric extension to induce passive adiabatic decompression melting of both asthenospheric and lithospheric upper mantle has been previously proposed to explain the origin of the anorogenic magmatism in the Mediterranean region (e.g., Lustrino and Wilson, 2007, and references therein). In this context, the 410 depression can be attributed to a compositional origin. Decompression melting in the overlying mantle preferentially extracts Fe, leaving a Mgenriched residue that should be transported into the mantle transition zone (Schmerr and Garnero, 2007). In a Mg-enriched mantle in the vicinity of the 410, the olivine-wadsleyite transition moves to higher pressures (Fei and Bertka, 1999) leading to the downward deflection of the 410 discontinuity beneath north Morocco and the northeast coast of Spain. The other group of models proposed to explain the anorogenic magmatism in the Mediterranean region requires active asthenospheric (or deeper) mantle convection, generally related to the Canary-Cape Verde plume (e.g., Lustrino and Wilson, 2007 and references therein). Plume-like models provide a thermal origin to the downward deflexion of the 410. The 410 depression (of about 15-20 km) can be attributed to a mantle temperature increase of about 161-215 K, which only affects the 410 discontinuity beneath north Morocco and the eastern coast of Spain. In this context, the Alboran slab would be an obstacle for the plume material, leading to an increase in the upper-mantle temperature towards the west of the Alboran slab. As a result the TZ becomes thinner beneath the Strait of Gibraltar, Gulf of Cadiz, and south of Portugal. We have no preferred model since the contribution of the Mg-enriched mantle to the depth of the 410 has the same sign as that of a higher temperature. Besides, as mentioned in Lustrino (2011), yet there are no evidences of the existence (or absence) of a physically continuous mantle source to support (or to question) the plume-like model. The lack of data in the southeast corner of our study region does not allow to physically relate the origin of the 410 downward deflection beneath the different areas.
The combined analysis of detailed upper mantle structure and high-resolution mapping of topography discontinuities over the whole Euro-Mediterranean region will help to further constrain the origin and nature of the anorogenic magmatism in the framework of the regional dynamics.

Conclusions
We have carried out a semi-automatic search for converted phases P410s and P660s beneath the Ibero-Maghrebian region using teleseismic data from 258 stations. The large data volume provided by the IberArray of the Spanish TopoIberia project has permitted us to considerably increase the resolution and to resolve new TZT and topography of discontinuities.
We have found lateral variations in the 410 and 660 topographies, which lead to TZT variations and correlate with known seismic anomalies in the Ibero-Maghrebian region. We have interpreted our results in terms of compositional and temperature variations inside the TZ. Our main conclusions are: • The thickened TZ beneath the Alboran Sea is consistent with the position of the high-velocity anomaly in the tomographic images related to the Betic-Alboran slab. Our results suggest that the slab is cold enough to induce downward deflection of the post-spinel transformation and uplift of the olivine-to-wadsleyite phase transition. • The anti-correlated depth of the 410 and 660 beneath the Alboran Sea have provided indirect evidences confirming the steep nature of the Betic-Alboran slab. • Significant lateral variation in the TZT beneath Granada coincides with the location of a nest of deep earthquakes. This variation has been interpreted as a high temperature gradient inside the TZ. • The presence of strong Pws phases before P410s and the depressed 660 in the Pyrenees and Balearic Sea are consistent with the presence of the Alpine-Tethys remnant slab. The cold stagnant slab can explain the 660 depression and the presence of a water concentration increase in the TZ due to dehydration of the stagnant slab can explain the strong reversed polarity signal before P410s. • The TZT beneath the western Moroccan region, the Gulf of Cadiz, and the southwest coast of Portugal may be controlled by the presence of multiple discontinuities at the base of the TZ. A compositional origin for the thickened TZ beneath the western Moroccan region has been proposed. The garnet-to-perovskite transition together with an Al-rich mantle can explain the downward deflection of the 660 beneath Morocco and the thicker TZ in this area. For the Gulf of Cadiz Fig. 11. Schematic explanation for the proposed models, where the topography of the 410 and 660 discontinuities is governed by the temperature inside the TZ. Unless indicated, the discontinuities are due to olivine phase transitions. and the southwest coast of Portugal a thermally thinned TZ has been proposed, which is consistent with the tomographic images from Monna et al. (2013). • The presence of melt atop the 410 may explain reversed polarity arrivals (Pws) from above this discontinuity. Plume material beneath the Gulf of Cadiz and subducted oceanic-lithosphere beneath the Alboran Sea and northeast Spain could provide the higher water content to explain the strong reversed polarity signals in these areas. • The spatial correlation between the active anorogenic magmatism and the 410 depression beneath Morocco and the northeast coast of Spain can be explained with both groups of models proposed for the origin of this magmatic activity. The estimated 410 depression in this area is consistent with a plume-like model with the stem of the plume away from the Mediterranean and affecting only the upper portion of the TZ, or with a Mg-enriched mantle due to decompression melting in the overlaying mantle, which shifts the 410 phase transition to higher pressures.