The Rif–Betics–Alboran region has been vital in the tectonic evolution of the western Mediterranean. Seismic images support the idea of continuous slab rollback being a prominent force in this region. However, the detailed slab structure and the physical mechanisms generating local deep (> 600 km) earthquakes remain unclear. Here, we analyze waveforms recorded from dense seismic networks above the deep earthquake beneath Granada in 2010 to study the slab structure. We discover a thin low‐velocity layer (LVL) at the base of the slab to explain both the long codas observed in Morocco and the secondary arrivals observed in Spain. This LVL indicates the presence of hydrous magnesium silicates extending to ∼600 km depth, which suggests that dehydration embrittlement promotes the occurrence of deep‐focus earthquakes. Our findings contradict the traditional slab model with the LVL sitting on the top of the slab, suggesting that the Alboran slab has been overturned.

Deep‐focus earthquakes with a depth of 300–700 km are mostly observed in subducted oceanic plates along convergent boundaries. At these depths with the extreme confining pressure and high temperatures, earthquakes are distinct from those at the shallower depths: significant nondouble‐couple mechanism, larger fracture energies, largely varied b‐values and aftershock sequences, and highly varied source durations and stress drops (see Zhan, 2020 for a review). Various physical mechanisms, including dehydration embrittlement, transformational faulting, and shear instabilities (Zhan, 2020 and references therein), have been proposed to explain deep earthquakes. However, seismological inquiries into the physical mechanism of deep earthquakes are impeded by the limited observations, in particular, the critical observations from local dense seismic arrays.

Among deep‐focus earthquakes occurred in different regions, those beneath southern Spain are unique. Since 1954, there have been five closely located large deep‐focus earthquakes with depths greater than 600 km beneath Granada (Fig. 1; Buforn et al., 2011). In contrast to deep‐focus earthquakes along the cold Tonga subduction that can have significant aftershocks (Wiens and Gilbert, 1996), the Spain deep‐focus earthquakes lack aftershocks, which is more similar to those beneath the Japan–Kuril region. Furthermore, no seismicity has been recorded between the depths of 150 and 600 km beneath southern Spain (Buforn et al., 2004). Although the origin of the Spain deep‐focus earthquakes is an open question, they are naturally linked to the subducted cold lithosphere. However, whether the high‐velocity structure is a subducted slab or delaminated continental lithosphere is still controversial, which adds another question about the evolution of the western Mediterranean.

Figure 1.

Map of seismic stations and the deep Granada earthquake (occurred on 11 April 2010, depth = 623 km, Mw=6.2) used in this study. (a) The stations are from the PICASSO seismic experiment (YB network) and the IberArray (IB network). The red triangles indicate the stations for studying the low‐velocity layer at the base of the slab. The blue inverted triangles in Spain indicate stations for studying anomalous P arrivals. The blue stars show the historic deep earthquakes between 1954 and 2010 (Buforn et al., 2011). (b) Two cross section AA′ and BB′ in panel (a) through the P‐wave tomography model of Bezada et al. (2013). The cyan dashed lines outline the hybrid slab model in Figure 2. The red and black lines show the ray paths for the P phase to the north and the S phase to the south, respectively.

Figure 1.

Map of seismic stations and the deep Granada earthquake (occurred on 11 April 2010, depth = 623 km, Mw=6.2) used in this study. (a) The stations are from the PICASSO seismic experiment (YB network) and the IberArray (IB network). The red triangles indicate the stations for studying the low‐velocity layer at the base of the slab. The blue inverted triangles in Spain indicate stations for studying anomalous P arrivals. The blue stars show the historic deep earthquakes between 1954 and 2010 (Buforn et al., 2011). (b) Two cross section AA′ and BB′ in panel (a) through the P‐wave tomography model of Bezada et al. (2013). The cyan dashed lines outline the hybrid slab model in Figure 2. The red and black lines show the ray paths for the P phase to the north and the S phase to the south, respectively.

The prevailing tectonic evolution models suggest that the western Mediterranean subduction system mainly retreated westward toward Gibraltar during the Miocene (Chertova et al., 2014). With the extensional collapse of the Alborán domain, the Betic–Rif arc forms its present dramatic arc geometry. However, which deep mantle processes, such as delamination, slab break‐off, convective removal of lithosphere, and slab rollback, are driving the extensional environment remain inconclusive (Platt et al., 2013). Although the recent detailed seismic tomographic models (Gutscher et al., 2002; Spakman and Wortel, 2004; Bezada et al., 2013; Palomeras et al., 2014; Villaseñor et al., 2015) confirm the existence of a subducted slab, the question of whether the slab belongs to the oceanic or continental crust lithosphere remains unanswered. Thus, distinguishing the origins of these deep mantle anomaly will be crucial for deciphering the dynamics of the slab–mantle interaction across the western Mediterranean.

Here, the densely spaced deployment of IberArray and PICASSO Array (Fig. 1) recorded the 11 April 2010, deep Spain earthquake (Mw 6.3, event depth of 623 km from Global Centroid Moment Tensor solution). The dense number of stations right above the deep‐focus earthquake provides an unprecedented opportunity to study both the physical mechanism of the earthquake and the detailed slab structure beneath the western Mediterranean. Through modeling the coda waves and extra P arrivals recorded by seismic stations in Morocco and Spain, respectively, we demonstrate the existence of hydrous magnesium silicates along the slab surface and surrounding the earthquake hypocenter, suggesting the important role of water in generating a deep‐focus earthquake. In addition, these data provide a new image of an overturned slab, which is important for better understanding the complicated subduction history in the western Mediterranean.

Seismograms for stations in Morocco (Fig. 1) show different coda waveforms depending on their locations. The southern stations located in the Atlas with epicentral distance Δ>3° show relatively simple waveforms (Fig. 2). In contrast, the stations near the Rif (red triangles in Fig. 1a) show strong coda waves after the direct P and S waves on all three components (Fig. 2, Fig. S1, available in the supplemental material to this article). On the SH component (Fig. 2), the coda last for ∼20 s. Such long coda arrivals are not present for teleseismic events (Fig. S2), suggesting that localized crustal structures beneath stations are not responsible for the strong coda. Thus, an anomalous mantle structure along the ray path is required for producing the long codas.

Figure 2.

Comparison between S data and synthetics for different models. The tangential data are selected from the stations indicated as the red triangles in Figure 1a. The red dashed line at 0 s indicates the predicted S arrival from the IASP91 model. From left to right, synthetics in panel (b) are calculated for the cross section AA′ (Fig. 1b) through the P‐wave tomography model of Bezada et al. (2013), a hybrid slab model, a hybrid slab model with a thin low‐velocity layer (LVL), a heterogeneous slab with randomly distributed scatterers, and a heterogeneous slab plus the LVL as displayed in panel (a). In the hybrid slab model, the δVS at the center of the slab is 4.5% and linearly decreases to 1.5% at the edge. In the models with the LVL, the LVL is at the base of the slab and has a width of 2 km and δVS of −10%. Note: For visualization purposes here, the thickness of the LVL is amplified from its true dimensions. The scatterers in the heterogeneous slab are characterized by a von Kármán distribution function with a down‐dip correlation distance a of 10 km and a perpendicular correlation distance b of 1 km. The scatterers have an elongated shape parallel to the slab dip. The standard deviation of the δVS (σVS) of the scatterers is 5%. Note the earthquake (red star) is located at the base of the slab as shown in the schematic in panel (a). Both data and synthetics are filtered to a frequency band of 0.5–5 Hz. The focal mechanism used in the simulations is from the Global Centroid Moment Tensor solutions (Global CMT; see Data and Resources) with strike, dip, and rake of 63°, 28°, and −39°, respectively.

Figure 2.

Comparison between S data and synthetics for different models. The tangential data are selected from the stations indicated as the red triangles in Figure 1a. The red dashed line at 0 s indicates the predicted S arrival from the IASP91 model. From left to right, synthetics in panel (b) are calculated for the cross section AA′ (Fig. 1b) through the P‐wave tomography model of Bezada et al. (2013), a hybrid slab model, a hybrid slab model with a thin low‐velocity layer (LVL), a heterogeneous slab with randomly distributed scatterers, and a heterogeneous slab plus the LVL as displayed in panel (a). In the hybrid slab model, the δVS at the center of the slab is 4.5% and linearly decreases to 1.5% at the edge. In the models with the LVL, the LVL is at the base of the slab and has a width of 2 km and δVS of −10%. Note: For visualization purposes here, the thickness of the LVL is amplified from its true dimensions. The scatterers in the heterogeneous slab are characterized by a von Kármán distribution function with a down‐dip correlation distance a of 10 km and a perpendicular correlation distance b of 1 km. The scatterers have an elongated shape parallel to the slab dip. The standard deviation of the δVS (σVS) of the scatterers is 5%. Note the earthquake (red star) is located at the base of the slab as shown in the schematic in panel (a). Both data and synthetics are filtered to a frequency band of 0.5–5 Hz. The focal mechanism used in the simulations is from the Global Centroid Moment Tensor solutions (Global CMT; see Data and Resources) with strike, dip, and rake of 63°, 28°, and −39°, respectively.

We first examine how accurately the tomographic model of Bezada et al. (2013) can explain the data at longer periods of 3–20 s, at which the overall properties of the slab can be represented. As displayed in Figure 1b, the southward ray path samples the slab edge and moves away from the slab with increasing epicentral distance. This is evident from the shift in the travel times of the direct S arrival (dTS), which gradually increases within the distance of 2°–2.5° and remains unchanged until to ∼3° (Fig. S3). Given the computing capability, we generate 2D synthetics along the cross section AA′ through the tomographic model (Fig. S4). We use a 2D finite‐difference scheme (Li et al., 2014) to generate the 2D synthetics accurate up to 8 Hz with the implementation of accurate point source mechanism.

Neither the original tomographic model nor the enhanced model (Fig. S3a,b), which multiplies the velocity perturbations in the original model by a factor of 3 (Fig. S3a), can explain the dTS. Furthermore, the smooth tomographic model has difficulty explaining the amplitude variation. In the data, the amplitude of SH drops at 2.3° distance and reaches the maximum at ∼3°, which is a strong indication of the multipathing effects due to sampling along the slab edge. Thus, we try hybrid slab models characterized by a simplified rectangular shape, encompassing both a vertical slab (Fig. S3a) and slabs featuring various dipping angles (Fig. S3b). Through tests for different source locations relative to the slab edges (Fig. S3c,d), we find that a southwest tilting of a subvertical slab structure with the earthquake located close to the southwest edge offers a reasonable fit of both the travel time and amplitude (Fig. S3d).

At higher frequencies of 0.5–5 Hz, the hybrid model produces a multipathed SH with two arrivals, but it does not produce strong coda arrivals (Fig. 2). In Figure S5, we further test different source locations and a vertical hybrid slab model, which proves inadequate for explaining the waveform complexity in the data. It has been previously found that elongated random heterogeneities inside the slab can produce high‐frequency codas efficiently. We use an identical random scatterer model as in the slabs beneath northern Japan and Calabria (Furumura and Kennett, 2005; Sun et al., 2014), and produce a strong coda for those stations that are located above the shallowest part of the continuous slab (Fig. 2). Because the event is located at the edge of the slab, the coda energy drops significantly when Δ>2.3° even with the scatterers present. In contrast, the strong coda can be seen up to 3° in the data. In Figures S6–S10, we test different scatterer models as well as the location of the event within the modelled slab. Notably, a vertical slab incorporating scatterers with a source located at the northeast edge of the slab (Fig. S11) produces coda similar to the data, but at a larger distance this model has difficulty to explain the travel times. Thus, we conclude that, with the event located close to the southwest edge of the slab, the scatterers in the slab can contribute to the development of the coda similar to the data, particularly at small distances. However, the station gap that is the result of the Alborán Sea (right above the slab) and lack of ocean bottom seismometers limits the resolution on the heterogeneities within the slab.

Transformational faulting has been the prevailing hypothesis for generating deep earthquakes (Green et al., 2010). Inside the cold slab, polymorphs of olivine can extend into the mantle transition zone forming a metastable olivine wedge (MOW), which hosts the occurrence of transformational faulting. The MOW as a low‐velocity structure has been identified in the Japan subduction zone (Iidaka and Suetsugu, 1992; Shen and Zhan, 2020). Here, we test slab models embedded with a MOW (Fig. S12). The low‐velocity MOW produces strong multipathed waveforms, which extend over large distances. At the distance range of 2°–3°, the later arrival sampling the core of the metastable olivine wedge becomes the dominant phase. Furthermore, no matter how we change the geometry and velocity perturbation of the MOW (Fig. S11), we cannot generate coda arrivals as in the data. Thus, a low‐velocity MOW here is not responsible for the strong coda.

A low‐velocity layer (LVL) atop the slab promotes the high frequency coda arrivals by trapping the waves in the LVL as a guided wave (Abers, 2000; Martin et al., 2003). We first test a model with an undulating low‐velocity serpentine layer on top of the subducting slab down to 600 km (Fig. S13), which is proposed to explain the ∼70 s coda at the frequency of 1 Hz observed for the deep earthquakes underneath Tonga (Savage, 2012). With a δVS of −10%, the undulating LVL models generate strong coda arrivals (Fig. S13). However, the initial SH arrival is much weaker than what is observed in the data, and the distance range showing strong coda extends to a much larger distance than is observed. A model with low‐velocity fault zones as in Garth and Rietbrock (2017) but extending to the mantle transition zone was also tested, which cannot predict the data either (Fig. S14).

We further test simplified layered LVL models defined by parameters of shear‐velocity perturbation (δVS), thickness (Th), and top depth from the surface (Du) (Fig. S4). An LVL model only extending into the shallow mantle (Du<100  km) can reproduce the coda behavior well. If the top of the LVL is present to a deeper depth of Du>200  km (Fig. S15), the guided wave energy leaks out from the LVL at those deeper depths and then spreads across a much broader distance range. For the same Th, the model with larger δVS generates stronger and longer codas (Fig. S16). A strong trade‐off exists between the δV and Th. For example, a model with Th of 0.5 km and δVS of −40% generates almost identical waveforms predicted from a model with Th of 2 km and δVS of −10% (Fig. S17). However, the shear‐velocity perturbation cannot be too small. For example, with δVS of −5%, an extremely thick LVL produces multipathed arrivals instead of the observed coda. Furthermore, the earthquake must be in or near the LVL to generate guided waves efficiently. If the earthquake is located outside the LVL by 5 km, the coda wave is much weaker (Fig. S18). We also test models by including both LVL and slab heterogeneities (Fig. 2), which better reproduces the coda at the distance of 2.5°–3°. However, such a model also produces strong coda at large distances. It is interesting to note that if we slightly shift the earthquake location into the interior of the slab (Fig. S19), the coda at distance larger than 3° is weaker, which is more representative of the data. By including both a LVL and scatterers, the high‐frequency codas of P and SV can be also explained (Fig. S20). It is worth noting that our primary objective here is to find a model efficiently captures the most significant features in the data instead of matching the data exactly, considering trade‐offs among different parameters. Nevertheless, a continuous thin LVL extending from the source region to the shallow mantle can explain the strong coda recorded on the stations in Morocco.

Stations in Spain recorded clear signals of two large subevents of the Granada deep earthquake on the direct P arrival (P1–2 in Fig. 3 and Fig. S21). Besides these two strong subevents, a strong extra arrival (P*) appears 0.5 s after the P1 at the stations north of the event near the nodal plane (red waveforms in Fig. 3). Because the P* is not coherently observed on most seismograms, both finite‐fault inversion from the teleseismic data and waveform analysis (Buforn et al., 2011; Bezada and Humphreys, 2012) do not identify the P* as a subevent. If the rupture propagates along a subhorizontal fault plane from east to west (Buforn et al., 2011), it is difficult to generate an azimuthally abrupt emergence of the P*.

Figure 3.

Anomalous P arrivals and the low‐velocity region surrounding the deep earthquake. (a) Displays the waveform records of the P waveforms on the vertical component. (b) Two representative waveforms as indicated by blue circles in panel (a). The data from station EMIN show two arrivals of the subevents P1 and P2. In contrast, station E071 has an extra arrival between the two P arrivals. The waveforms in panel (a) showing apparent extra arrival are marked with red, which is focused in an azimuthal range of ∼−30° to 30°. (c) The extra arrival is modeled as an S‐to‐P converted phase at the interface between the LVL (orange color) and the slab (blue color). (d) Comparison between the stacked P data (black traces) in the azimuth of 0°–10° and synthetics (red traces) for different LVL models. Here, Th describes the thickness of the LVL near the earthquake, and D is the lateral distance from the earthquake to the interface between the LVL and the slab as shown in panel (c). The number at the top of the traces represents the cross‐correlation coefficient (CC) between the data and synthetics. The LVL model with the highest CC has a thickness of 6 km and a D of 4.8 km. In the waveform modeling, we use different mechanisms for the two subevents, which correspond to E2 and E3 in Bezada and Humphreys (2012) with strike/dip/rake of 64°/24°/−45° and 79°/27°/−23°, respectively.

Figure 3.

Anomalous P arrivals and the low‐velocity region surrounding the deep earthquake. (a) Displays the waveform records of the P waveforms on the vertical component. (b) Two representative waveforms as indicated by blue circles in panel (a). The data from station EMIN show two arrivals of the subevents P1 and P2. In contrast, station E071 has an extra arrival between the two P arrivals. The waveforms in panel (a) showing apparent extra arrival are marked with red, which is focused in an azimuthal range of ∼−30° to 30°. (c) The extra arrival is modeled as an S‐to‐P converted phase at the interface between the LVL (orange color) and the slab (blue color). (d) Comparison between the stacked P data (black traces) in the azimuth of 0°–10° and synthetics (red traces) for different LVL models. Here, Th describes the thickness of the LVL near the earthquake, and D is the lateral distance from the earthquake to the interface between the LVL and the slab as shown in panel (c). The number at the top of the traces represents the cross‐correlation coefficient (CC) between the data and synthetics. The LVL model with the highest CC has a thickness of 6 km and a D of 4.8 km. In the waveform modeling, we use different mechanisms for the two subevents, which correspond to E2 and E3 in Bezada and Humphreys (2012) with strike/dip/rake of 64°/24°/−45° and 79°/27°/−23°, respectively.

The P* is also noted in Bezada and Humphreys (2012), which they interpreted as wave interactions with the 660 km discontinuity, the upper‐mantle structure, and the shallow structure beneath the stations. We test different properties of the 660 km discontinuity and fail to produce the strong P* as a reflection from that discontinuity. Along the northward ray path (Fig. 1b), no apparent strong, small‐scale upper‐mantle anomalies are presented in the tomographic models, which can produce significant multipathing effects. Furthermore, the geographical distribution of stations showing strong P* in Figure 3a does not correlate with any major tectonic structures in the Spanish continental lithosphere. We also examine both regional and teleseismic events, and find no extra arrivals recorded at these stations, which suggests that P* is not caused by shallow crustal structures beneath the stations but has a deeper origin near the source.

We suggest that a locally converted SP wave at the source region produces the P* phase (Matsuzawa et al., 1990). The sharp velocity contrast across the interface between the LVL and slab allows the shear wave converting to P waves and forming the P* as a later arrival (Fig. 3c), which can also explain the prominent P* presented near the P‐wave nodal plane. It is interesting to note that we also observe similar P* for a deep event that occurred in the northwestern Pacific region (Fig. S22), which also appears strongly near the P‐wave nodal. Although we cannot totally exclude the P* as a subevent, the SP wave is probable, considering the existence of the LVL at the slab surface.

Despite the uncertainties of the derived focal mechanisms of subevents P1–2 and the tilting angle of the interface between the slab and the LVL, we group the data in 20° azimuthal bins (Fig. 3c, Fig. S23) and model the stacked P data by assuming that the LVL interface is parallel to the slab. As displayed in Figure 3d, a larger δVP produces a stronger P*, and a larger distance (D) between the event and the interface generates a delayed P*. Thus, the amplitude and travel time of the P* provide good constraints on both the velocity and size of the LVL surrounding the source. By calculating the cross‐correlation coefficients between the data and synthetics for different δVP and D, we conclude that a localized LVL with δVP of 15%–20% and D of 4–6 km can best predict both amplitude and travel time of the P* for the three azimuthal bins (Fig. 3d, Fig. S23).

Implication of the origin of the LVL

By modeling the strong coda at the stations in Morocco and the P* at the stations in Spain, we confirm the existence of the LVL along the top of slab surface as well as surrounding the Spain deep‐focus earthquake. At depths shallower than 200 km, a seismic low‐velocity layer on the top of subducting slabs has been routinely observed globally with guided waves, which has been interpreted as partially metamorphosed subducted crust (Abers, 2000). However, the LVL imaged here extends continuously to the mantle transition zone, suggesting a different mechanism. Li et al. (2018) suggest that laminated fabric or aligned melt pockets parallel to the slab interface can produce a shear‐wave anisotropy of ∼25%, which can explain the continuous LVL we have imaged and its velocity reduction up to ∼20%. However, we also see strong codas following both P and SV (Fig. S1) for stations in Morocco, which indicates that the LVL has a strong velocity reduction for P and SV, and presents a weak anisotropy.

Alternatively, the existence of hydrated phases from the mantle lithosphere, that is, serpentine and phase A may provide the most feasible explanation of the LVL. For the LVL attached to the slab surface, 62% mantle serpentinization occurs if a large VS reduction of 40% corresponding to a δVP of −26% (assuming δVs/δVP=1.5) is used, and the water content in serpentine is 8.1 wt% with w(%)0.31δVP (Carlson and Miller, 2003). This inferred mantle serpentinization ratio, assuming the water content of the serpentine is ∼13% (Carlson and Miller, 2003), is similar to the 60% calculated by Savage (2012). In contrast, using a low‐VS reduction of 10% and VP of −6.7% generates 16% mantle serpentinization, which is similar to the value of 15%–30% of the shallow mantle from seismic inversions (van Avendonk et al., 2011).

For the LVL surrounding the source, we have modeled a VP drop of 15%–20%. Such a large velocity drop is difficult to be explained by metastable olivine, which has a VP drop smaller than 10% (Mao et al., 2015). Thus, it is possible that the localized LVL is due to hydrous magnesium silicates, which suggests significant water surrounding the earthquake hypocenter. To retain such a high‐water content all the way to the mantle transition zone, the subducted plate should be cold to have stable hydrous silicates from the depth of ∼300 km to the mantle transition zone (Komabayashi et al., 2004). Thus, with a relatively young sea floor age in the western Mediterranean, the subduction here is at least at a moderate speed, similar to ∼70 mm/yr in northwestern Pacific subduction, to remain cold (van Keken et al., 2011).

Water and deep earthquake

The prevalence of using dehydration embrittlement to explain the mechanism for deep earthquakes is largely impeded by the insufficient water transported below 300 km through the main body of slab (Green et al., 2010). However, the feasibility of a thin layer atop the slab transporting water into greater depths remains unclear (Shen and Zhan, 2020). The observations of the LVL here and in the subducted lithosphere beneath Tonga–Fiji suggest that significant water can be carried to the mantle transition zone via subduction. Thus, we suggest that dehydration embrittlement can play a pivotal role in the occurrence of deep‐focus earthquakes beneath Spain. In addition, the transfer of stress or released fluid to the brittle surrounding regions may also trigger seismicity (Ferrand et al., 2017).

The dehydration cause for the deep event can also explain the distribution of other deep earthquakes beneath Spain, the scarce but clustered deep‐focus earthquakes, and the absence of intermediate‐depth seismicity. The deep‐focus earthquakes beneath the Peru–Brazil border (Zahradník et al., 2017) and northern Chile (Omori et al., 2004) show similar patterns. Thus, to produce the deep seismicity at ∼600 km, a cold slab with a steep geotherm is required to produce dehydration, possibly involving hydrous magnesium silicates like phase D to superhydrous phase B (shB) or series of dehydration reactions from shB (Komabayashi et al., 2004; Omori et al., 2004). However, the precise dehydration reactions in the mantle transition zone and the conditions for their occurrences remain elusive. Such a steep cold geotherm is also supported by the imaged subvertical slab and the continuous LVL. However, the subducted slab beneath the western Mediterranean cannot be as cold as the Tonga slab, which could have dehydration over a wider pressure range and has strong seismicity across many depths. Furthermore, the existence of water would enhance the transformation kinetics of olivine (Hosoya et al., 2005), which is also detrimental to the presence of MOW, but Ishii and Ohtani (2021) show that a MOW can exist for a wet slab. Thus, we cannot completely rule out transformational faulting nor shear instabilities, as the possible mechanism of deep‐focus earthquakes. In any case, hydrous magnesium silicates and associated dehydration may play a critical role in the generation of deep earthquakes, and the region hosting the clustered deep‐focus earthquakes beneath Spain may represent a “water tank” with limited size.

An overturned slab beneath the Betics

The existence of the continuous LVL confirms that the slab beneath the Betics of southern Spain is a subducted oceanic lithosphere. Thus, the formation of the Beltic–Rif arc must be due to subduction processes, either by slab rollback or slab break‐off. Unlike normal subduction with a low‐velocity layer on top of the slab surface, this slab (Fig. 4) has the low‐velocity layer “underneath” the slab surface dipping to the northeast. To produce such an overturned slab geometry, one may expect extreme scenarios such as northward subduction with significant trench advance. Although a northward subduction beneath the Betics is presented in the reconstruction model by Faccenna et al. (2004) with an initial subduction zone extended from the Baleares margin to Gibraltar in the Early Oligocene, such a model is less favored by 3D slab‐evolution modeling results (Chertova et al., 2014). Furthermore, despite of different subduction geometries, instead of trench advance, subduction rollback is favored by most tectonic evolution models of the Mediterranean. Alternatively, we propose that the current geometry of the slab beneath Betics is an overturned slab (Fig. 4), which initially subducts northward, and then is formed subsequently with slab rollback and docking on both Iberian and African continental margins. As shown in Moresi et al. (2014), subduction involving both rollback and continental accretion can generate complicated 3D structures, including overturned slabs. Nonetheless, a refined slab structure will be essential for better constraining the tectonic evolution history of the Mediterranean subduction systems.

Figure 4.

Schematic illustration of the slab evolution in the western Mediterranean region. (a) The retreat of the trench, indicated by the red solid line, due to the slab rollback. The magenta and cyan dashed lines depict the slab surface at the depth of 200 and 600 km, respectively. (b) The slab morphology along the cross‐section A. Note that the present slab is overturned, which causes the slab’s upper surface and LVL appearing as the bottom of the high‐velocity slab.

Figure 4.

Schematic illustration of the slab evolution in the western Mediterranean region. (a) The retreat of the trench, indicated by the red solid line, due to the slab rollback. The magenta and cyan dashed lines depict the slab surface at the depth of 200 and 600 km, respectively. (b) The slab morphology along the cross‐section A. Note that the present slab is overturned, which causes the slab’s upper surface and LVL appearing as the bottom of the high‐velocity slab.

All seismic data are archived and openly available at the Incorporated Research Institutions for Seismology Data Management Center (IRIS‐DMC; www.iris.edu/ds) and the Observatories and Research Facilities for European Seismology (ORFEUS) Data Center (www.orfeus-eu.org/fdsnws/dataselect/1/), under station code XB for the PICASSO Array (doi: 10.7914/SN/XB_2009), and IB for the IberArray (doi: 10.7914/SN/IB). All maps in this article were produced using Generic Mapping Tools (GMT; www.soest.hawaii.edu/gmt) developed by Paul Wessel and Walter H. F. The other relevant data to this article were available at https://www.globalcmt.org. All websites were last accessed in June 2023. The supplemental material comprises Figures S1 to S22: Figures S1 and S2 display PSV data of the deep earthquake studied and records of teleseismic events; Figures S3–S19 present synthetic results for various slab and low velocity layer models, along with variations in event location; Figures S20 and S21 show detailed waveforms for deep earthquakes beneath Spain and northeastern China, respectively; and Figure S22 illustrates the detailed modeling of the SP phase.

The authors acknowledge that they are no conflicts of interest recorded.

The authors thank Editor‐in‐Chief Keith D. Koper, Associate Editor Pascal Audet, and two anonymous reviewers for the constructive comments that help significantly improve the article. The authors thank the Incorporated Research Institutions for Seismology Data Management Center (IRIS) Data Center and Observatories and Research Facilities for European Seismology (ORFEUS) for providing the seismic data. All maps in this article were produced using Generic Mapping Tools (GMT) developed by Paul Wessel and H. F. Walter. The authors appreciate the Supercomputing Center of USTC for high‐performance computing services. Daoyuan Sun is supported by National Natural Science Foundation of China, Grant Number 41722401, B‐type Strategic Priority Program of the Chinese Academy of Sciences, Grant Numbers XDB41000000 and XDB42000000, and the Fundamental Research Funds for the Central Universities (Grant Number WK2080000144). The data for network XB were collected via funding by National Science Foundation Continental Dynamics Program under Grant Number EAR‐0809023 and Meghan S. Miller and Daoyuan Sun were supported by the National Science Foundation Geophysics Program under Grant Number EAR‐1345015 and National Science Foundation CAREER Award Number EAR‐1054638.

1.
Abers
G. A.
2000
.
Hydrated subducted crust at 100‐250 km depth
,
Earth Planet. Sci. Lett.
  176,
323
330
.
2.
Bezada
M. J.
, and
Humphreys
E. D.
2012
.
Contrasting rupture processes during the April 11, 2010 deep‐focus earthquake beneath Granada, Spain
,
Earth Planet. Sci. Lett.
  353,
38
46
.
3.
Bezada
M. J.
Humphreys
E. D.
Toomey
D. R.
Harnafi
M.
Davila
J. M.
, and
Gallart
J.
2013
.
Evidence for slab rollback in westernmost Mediterranean from improved upper mantle imaging
,
Earth Planet. Sci. Lett.
  368,
51
60
.
4.
Buforn
E.
Bezzeghoud
M.
UdÍas
A.
, and
Pro
C.
2004
.
Seismic sources on the Iberia‐African plate boundary and their tectonic implications
,
Pure Appl. Geophys.
  161,
623
646
.
5.
Buforn
E.
Pro
C.
Cesca
S.
Udias
A.
, and
del Fresno
C.
2011
.
The 2010 Granada, Spain, deep earthquake
,
Bull. Seismol. Soc. Am.
  101,
2418
2430
.
6.
Carlson
R. L.
, and
Miller
D. J.
2003
.
Mantle wedge water contents estimated from seismic velocities in partially serpentinized peridotites
,
Geophys. Res. Lett.
  30, 1250, doi: .
7.
Chertova
M. V.
Spakman
W.
Geenen
T.
van den Berg
A. P.
, and
van Hinsbergen
D. J. J.
2014
.
Underpinning tectonic reconstructions of the western Mediterranean region with dynamic slab evolution from 3‐D numerical modeling: Western Mediterranean slab evolution
,
J. Geophys. Res.
  119,
5876
5902
.
8.
Faccenna
C.
Piromallo
C.
Crespo‐Blanc
A.
Jolivet
L.
, and
Rossetti
F.
2004
.
Lateral slab deformation and the origin of the western Mediterranean arcs
,
Tectonics
  23, TC1012, doi: .
9.
Ferrand
T. P.
Hilairet
N.
Incel
S.
Deldicque
D.
Labrousse
L.
Gasc
J.
Renner
J.
Wang
Y.
Green
H. W.
II
, and
Schubnel
A.
2017
.
Dehydration‐driven stress transfer triggers intermediate‐depth earthquakes
,
Nat. Commun.
  8, 15247, doi: .
10.
Furumura
T.
, and
Kennett
B. L. N.
2005
.
Subduction zone guided waves and the heterogeneity structure of the subducted plate: Intensity anomalies in northern Japan
,
J. Geophys. Res.
  110, no. B10, doi: .
11.
Garth
T.
, and
Rietbrock
A.
2017
.
Constraining the hydration of the subducting Nazca plate beneath northern Chile using subduction zone guided waves
,
Earth Planet. Sci. Lett.
  474,
237
247
.
12.
Green
H. W.
II
Chen
W.‐P.
, and
Brudzinski
M. R.
2010
.
Seismic evidence of negligible water carried below 400‐km depth in subducting lithosphere
,
Nature
  467,
828
831
.
13.
Gutscher
M. A.
Malod
J.
Rehault
J. P.
Contrucci
I.
Klingelhoefer
F.
Mendes‐Victor
L.
, and
Spakman
W.
2002
.
Evidence for active subduction beneath Gibraltar
,
Geology
  30, 1071, doi: .
14.
Hosoya
T.
Kubo
T.
Ohtani
E.
Sano
A.
, and
Funakoshi
K.‐I.
2005
.
Water controls the fields of metastable olivine in cold subducting slabs
,
Geophys. Res. Lett.
  32, L17305, doi: .
15.
Iidaka
T.
, and
Suetsugu
D.
1992
.
Seismological evidence for metastable olivine inside a subducting slab
,
Nature
  356,
593
595
.
16.
Ishii
T.
, and
Ohtani
E.
2021
.
Dry metastable olivine and slab deformation in a wet subducting slab
,
Nature Geosci.
  14,
526
530
.
17.
Komabayashi
T.
Omori
S.
, and
Maruyama
S.
2004
.
Petrogenetic grid in the system MgO‐SiO2‐H2O up to 30 GPa, 1600°C: Applications to hydrous peridotite subducting into the Earth’s deep interior
,
J. Geophys. Res.
  109, no. 
B3
, doi: .
18.
Li
D.
Helmberger
D.
Clayton
R. W.
, and
Sun
D.
2014
.
Global synthetic seismograms using a 2‐D finite‐difference method
,
Geophys. J. Int.
  197,
1166
1183
, doi: .
19.
Li
J.
Zheng
Y.
Thomsen
L.
Lapen
T. J.
, and
Fang
X.
2018
.
Deep earthquakes in subducting slabs hosted in highly anisotropic rock fabric
,
Nature Geosci.
  11,
696
700
.
20.
Mao
Z.
Fan
D.
Lin
J.‐F.
Yang
J.
Tkachev
S. N.
Zhuravlev
K.
, and
Prakapenka
V. B.
2015
.
Elasticity of single‐crystal olivine at high pressures and temperatures
,
Earth Planet. Sci. Lett.
  426,
204
215
.
21.
Martin
S.
Rietbrock
A.
Haberland
C.
, and
Asch
G.
2003
.
Guided waves propagating in subducted oceanic crust
,
J. Geophys. Res.
  108, no. 
B11
, 2536, doi: .
22.
Matsuzawa
T.
Kono
T.
Hasegawa
A.
, and
Takagi
A.
1990
.
Subducting plate boundary beneath the northeastern Japan arc estimated from SP converted waves
,
Tectonophysics
  181,
123
133
.
23.
Moresi
L.
Betts
P. G.
Miller
M. S.
, and
Cayley
R. A.
2014
.
Dynamics of continental accretion
,
Nature
  508,
245
248
.
24.
Omori
S.
Komabayashi
T.
, and
Maruyama
S.
2004
.
Dehydration and earthquakes in the subducting slab: empirical link in intermediate and deep seismic zones
,
Phys. Earth Planet. In.
  146,
297
311
.
25.
Palomeras
I.
Thurner
S.
Levander
A.
Liu
K.
Villasenor
A.
Carbonell
R.
, and
Harnafi
M.
2014
.
Finite‐frequency Rayleigh wave tomography of the western Mediterranean: Mapping its lithospheric structure
,
Geochem. Geophys. Geosys.
  15,
140
160
.
26.
Platt
J. P.
Behr
W. M.
Johanesen
K.
, and
Williams
J. R.
2013
.
The Betic‐Rif Arc and its orogenic hinterland: A review
,
Annu. Rev. Earth Planet. Sci.
  41,
313
357
.
27.
Savage
B.
2012
.
Seismic constraints on the water flux delivered to the deep Earth by subduction
,
Geology
  40,
235
238
.
28.
Shen
Z.
, and
Zhan
Z.
2020
.
Metastable olivine wedge beneath the Japan Sea imaged by seismic interferometry
,
Geophys. Res. Lett.
  47, no. 
6
, doi: .
29.
Spakman
W.
, and
Wortel
R.
2004
.
A Tomographic View on Western Mediterranean Geodynamics
 ,
Springer Berlin Heidelberg
,
Berlin, Heidelberg
.
30.
Sun
D.
Miller
M. S.
Piana Agostinetti
N.
Asimow
P. D.
, and
Li
D.
2014
.
High frequency seismic waves and slab structures beneath Italy
,
Earth Planet. Sci. Lett.
  391,
212
223
.
31.
van Avendonk
H. J. A.
Holbrook
W. S.
Lizarralde
D.
, and
Denyer
P.
2011
.
Structure and serpentinization of the subducting Cocos plate offshore Nicaragua and Costa Rica
,
Geochem. Geophys. Geosys.
  12, Q06009, doi: .
32.
van Keken
P. E.
Hacker
B. R.
Syracuse
E. M.
, and
Abers
G. A.
2011
.
Subduction factory: 4. Depth‐dependent flux of H2O from subducting slabs worldwide
,
J. Geophys. Res.
  116, no. 
B14
, doi: .
33.
Villaseñor
A.
Chevrot
S.
Harnafi
M.
Gallart
J.
Pazos
A.
Serrano
I.
Córdoba
D.
Pulgar
J. A.
, and
Ibarra
P.
2015
.
Subduction and volcanism in the Iberia–North Africa collision zone from tomographic images of the upper mantle
,
Tectonophysics
  663,
238
249
.
34.
Wiens
D. A.
, and
Gilbert
H. J.
1996
.
Effect of slab temperature on deep‐earthquake aftershock productivity and magnitude‐frequency relations
,
Nature
  384,
153
156
.
35.
Zahradník
J.
Čížková
H.
Bina
C. R.
Sokos
E.
Janský
J.
Tavera
H.
, and
Carvalho
J.
2017
.
A recent deep earthquake doublet in light of long‐term evolution of Nazca subduction
,
Sci. Rep.
  7, 45153, doi: .
36.
Zhan
Z.
2020
.
Mechanisms and implications of deep earthquakes
,
Annu. Rev. Earth Planet. Sci.
  48,
147
174
.
This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Supplementary data