Revealing the Secrets of the Western Mediterranean: A Deep Earthquake and the Overturned Slab Open Access

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.

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.

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<100km$⁠) can reproduce the coda behavior well. If the top of the LVL is present to a deeper depth of $Du>200km$ (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*.

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

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

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.

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.

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 P‐SV 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.