Deep-focus earthquakes occur in subducting slabs of lithosphere at depths greater than 300 km. The mechanisms of these earthquakes must be fundamentally different from the brittle failure associated with earthquakes on the shallow (<70 km) megathrust plate interface because of the vast differences in pressure, temperature, and composition (Frohlich, 2006; Green and Houston, 1995; Kirby et al., 1996). Three previously proposed rupture mechanisms for deep-focus earthquakes include a phase transformation in the olivine polymorphs (Green and Burnley, 1989; Kirby et al., 1991; Green et al., 2015), thermal shear instabilities (Kanamori et al., 1998; Wiens, 2001), and embrittlement resulting from devolatilization (e.g., dehydration) of volatile-bearing phases present over a suitable spatial extent (Jung et al., 2004). The last mechanism is widely invoked to explain intermediate-depth earthquakes between 70 and 300 km (Peacock, 2001; Hacker et al., 2003).

The May 24, 2013, Mw 8.3, 609 km deep, Sea of Okhotsk earthquake provides a new perspective in understanding the rupture mechanism of the deep-focus earthquakes on account of its unusually large magnitude and the rapidly increased density of seismic stations at the Earth’s surface. Previous studies of the rupture process of this earthquake revealed important aspects of the earthquake’s mechanics, but did not agree with each other in all important respects (Wei et al., 2013; Ye et al., 2013; Chen et al., 2014; Meng et al., 2014; Zhan et al., 2014). As found by Ruiz et al. (2017), it is of importance to identify the correct orientation of the causative faults of the deep earthquakes since incorrect assumption of the faults’ orientation could result in biased parameters, and thus may lead to misinterpretation to the mechanisms of the deep earthquakes. The 2013 Okhotsk earthquake was found to rupture sequentially on a subvertical plane by inverting teleseismic P, SH, and depth phases (pP and sSH) (Chen et al., 2014), but was discovered to rupture a subhorizontal fault plane using teleseismic P waves (Wei et al., 2013; Ye et al., 2013). Moreover, the subhorizontal fault as the causative fault is also supported by the P and pP-wave single-array back-projection (BP) imaging (Meng et al., 2014). In addition, the single-array back projection typically has the uncertainties of the subevents approximately 20 km (Zhang et al., 2015). In this study, the multi-array BP imaging is used to dramatically increase the resolution of the images as detailed in A Synthetic Test for the Causative Fault’s Orientation and BP Images’ Resolution.

To infer a comprehensive model for the rupture of the 2013 deep Okhotsk earthquake, we combined two different rupture-assessment methods: one is the multi-array BP method (Kiser and Ishii, 2012; Zhang et al., 2016), modified from the BP method (Ishii et al., 2005; Xu et al., 2009; Zhang and Ge, 2010; Zhang et al., 2011), and the other one is the multi-subevent inversion (Kikuchi and Kanamori, 1991). The back-projection analyses were carried out in two different frequency bands (Koper et al., 2011; Wang and Mori, 2011; Lay et al., 2012; Yao et al., 2013; Fan and Shearer, 2015) using data from USArray (TA) and European seismic network (EU) (Figure 1), while the multi-subevent inversion was performed using the data from the Global Seismic network (GSN) (Figure 2). Next, we integrated the back-projection and multi-subevent models with the local curved strike of the subducting Pacific slab. Additionally, to identify the orientation of the causative fault of the earthquake and resolution of the BP images, a synthetic test composed of two seismic point sources has been performed. Finally, we discuss the source properties of the earthquake and their implications for the mechanics of the earthquake on the basis of the resulting rupture models.

To determine the orientation of the causative fault (subhorizontal or subvertical), a BP synthetic test was carried out for the arrays TA and EU. The P waves were synthesized from two seismic point sources separating 20 km vertically with a 5-s bursting difference using the orthonormal propagator method (Wang, 1999) and were filtered in the frequency band of 0.1–0.5 Hz as shown in Figure 6. We performed the single-array BP as above mentioned for both arrays and then combined them. The BP images of the two point sources were derived as shown in Figure 7. The second source is imaged right at the same location of the first one in the combined BP image. This indicates that the vertical resolvability of the multi-array BP imaging is too low to discriminate the depth differences of the subevents, but the depth differences of the subevents have no effect on the multi-array BP imaging on the horizontal plane. If the rupture of the earthquake is on the subvertical plane, the width of the causative plane could be much greater than 121 km, given the dip angle of approximately 80^o from GCMT and the rupture width of ∼21 km projected to the horizontal plane during the second rupture stage (Figure 4). The shallow intraslab M > 8 earthquakes such as 1933 Mw 8.5 Sanriku earthquake (Okal et al., 2016) and 2017 Mw 8.2 Chiapas, Mexco earthquake (Zhang and Brudzinski, 2019) typically have rupture widths less than 100 km. Moreover, compared to the shallow earthquakes, the deep-focus earthquakes have faster rise times (Houston and Williams, 1991) and hence compacter rupture dimensions. The extremely wide rupture seems to be unlikely for the deep-focus earthquake. Therefore, the orientation of the causative fault of the 2013 Okhotsk earthquake is likely to be subhorizontal.

Besides confining the orientation of the causative fault, the synthetic test shows the resolution of the BP images. For the TA and EU, the second source is imaged ∼20 km northeast and ∼5 km south to the real location, respectively. However, for the combined image, the second source is imaged right to its real location. This indicates that the multi-array BP images have a high resolution less than 1 km, which is much better than those from the single-array BP. This is the reason why we take multi-array, multi-band, BP imaging in this study.

Our multi-array, multi-frequency-band back-projected rupture imaging and the inverted multiple subevents show that the deep earthquake’s asymmetrical, bilateral rupture occurred within 30 s on a subhorizontal, curved fault plane (Figure 5A). The rupture could be separated into three principal segments at rupture velocities of 2.7–4.6 km/s, 47–84% of shear velocity (Fukao and Obayashi, 2013) at the focal depth (Figures 4, 5B). The rupture velocity is constrained from the multi-array, multi-frequency BP. The first segment ruptured northwards and somewhat down-dip on the subhorizontal fault plane at a high rupture velocity of 4.3 km/s, 78% of shear velocity at the focal depth. Unsurprisingly, this segment radiated strong high-frequency seismic energy. After a transitional second segment ruptured, the third segment started on the other (south) side of the nucleation point and ruptured in the reverse direction: southwards and somewhat up-dip at a velocity of 2.7 km/s, 47% of shear velocity at the focal depth, while radiating low-frequency seismic energy. For the first and third segments, the considerable slab-normal components of the rupture are 23 and 31 km wide, which are estimated from the multi-subevent inversion. The slab-normal component of the second segment is constrained to 19 km using the rupture length of the BP image. The wide slab-normal rupture has also been observed by Meng et al. (2014). Each of the rupture distances is wider than the thickness of the metastable olivine wedge (MOW) (less than 15 km estimated by Kirby et al. (1996) or approximately 15 km constrained by Marone and Liu (1997)) using numerical modeling unless the slab at the focal depth is contorted. Additionally, if there exists the MOW, the nucleation of the earthquake would occur near the bottom of the MOW as shown in the inset of Figure 8 since there stores more strain energy for the occurrence of large earthquakes. Besides the high resolution of the multi-array, multi-frequency BP images, the location of the MOW relative to the rupture of the earthquake has an uncertainty less than 5 km.

To confirm the rupture beyond the possible metastable olivine wedge, a source model with 12 subevents (Supplementary Figure 5B) is inverted on a curved plane with the method of Kikuchi and Kanamori (1991) using the data in Figure 2B. The curved plane is designed parallel to the strike of the subducting slab. The resulting misfit of the waveforms is 0.56, which is much larger than that derived from the final multi-subevent model in Figure 5A and Supplementary Figure 5A. The better fitting of waveforms for the final multi-subevent model can be seen in Supplementary Figure 4, 6. This emphasizes the necessity of the rupture beyond the possible MOW.

Intriguingly, the rupture in the first 13 s extends northwards of the hypocenter, followed by the fourth subevent that bursts 44 km southwest of the epicenter with an azimuth of approximately 200^o, without subevents to spatially connect the two energy bursts (Figure 8A). The body wave incident on the location of the 13-secend subevent has a takeoff angle of approximately 90^o. Given the subhorizontal fault (strike 189^o, dip 11^o, and rake −93^o), the P-to-S amplitude ratio is only 0.082 according to radiation pattern (Aki and Richards, 2002), so triggering of the rupture by p waves is unlikely. The shear velocity at the depth of the hypocenter (609 km) is approximately 5.5 km/s in the velocity model AK135 (Kennett et al., 1995). The first S arrival at the 13-s subevent from the hypocenter is 8 s after the commence of the earthquake. The onset of the post-13-s subevent to the southwest is consistent with the arrival of S-waves from the preceding subevents. The appropriate timing and large amplitude of the S-waves thus indicate that the occurrence of the post-13-s subevent is dynamically triggered by the running S-waves. The coseismic dynamic triggering was also discovered in the 2012 Mw 8.6 Sumatra offshore earthquake (Zhang et al., 2012).

The rupture models in this study are basically consistent with those using single-array back-projection methods and/or inversion techniques (Wei et al., 2013; Ye et al., 2013; Meng et al., 2014; Zhan et al., 2014). For the fault plane, Wei et al. (2013) exclude the near-vertical fault plane based on the strong directivity effect along the azimuth of 165^o, which significantly deviates from the strike (12^o given by the GCMT) of this fault plane. The trending of subevents significantly off the strike of the subvertical fault plane leads Zhan et al. (2014) to choose the subhorizontal fault as the causative fault. Also, the subhorizontal fault plane is preferred by Ye et al. (2013) according to the single-array back-projection models and finite fault models. Additionally, the bilateral rupture in our rupture models agrees with another back-projection study by Meng et al. (2014).

The second largest subevent ruptures backward to the north relative to the largest subevent (Figure 5A). The rupture feature is confirmed by removing the two subevents in the new inversion model, respectively. The synthetic waveforms are significantly affected by the two subevents and show that the largest subevent radiates the first amplitude peak (Supplementary Figure 7), and the second largest one generates the second amplitude peak (Supplementary Figure 8). The locations and rupture times of the largest two energy-releasing peaks in our results agree reasonably well with other groups’ such as Zhan et al. (2014). This feature may be explained as a barrier, corresponding to the second largest subevent, is unbroken when the rupture first reaches it but eventually breaks down due to stress transferring. This intriguing rupture reversal has also been observed by prior studies (Wei et al., 2013; Zhan et al., 2014).

A profound difference between the models obtained in this study and by Ye et al. (2013) and Wei et al. (2013) is that the southernmost rupture is spatially shorter, resulting in this study’s more compact rupture area and lower rupture velocity. The compact rupture area in our model is the result of jointly back projecting seismic data from two arrays at different azimuths, reducing artifacts related to “smearing” of energy along the event-array azimuth (Zhang et al., 2016; Zhang et al., 2017; Zhang and Ge, 2017; Zhang and Brudzinski, 2019). The source model (Supplementary Figure 5C) by Zhan et al. (2014) lacks subevents in the northernmost, down-dip portion of our rupture models probably due to inversion of the displacement waveforms for a limited number of subevents.

The estimated extent of Okhotsk rupture does not favor simple transformational faulting in metastable olivine as the primary instability mechanism, because the rupture extent exceeds the likely plausible thermo-spatial extent of metastable olivine in the slab. Likewise, both the large seismic efficiency of 0.6 determined by Ye et al. (2013) and fast rupture velocities in the first two rupture segments preclude (Tibi et al., 2003) a thermal shear instability (Ogawa, 1987; Hobbs and Ord, 1988) as the rupture mechanism. The simplest explanation of the earthquake would be that dehydration embrittlement, associated with intermediate-depth earthquakes, still operates at these great depths (Omori et al., 2004). Three more speculative seismogenic options include: 1) another form of transformational faulting, perhaps transformation of metastable pyroxene to metastable akimotoite; 2) nucleation by transformational faulting in the metastable olivine wedge that propagates into surrounding rock by shear melting (Green, 2017), perhaps aided by latent heat release (Bina, 1998); 3) unusual slab contortions that completely enclose the rupture plane within a flattened MOW (Stein, 1995).

Option 1) may seem plausible given its exothermicity (Green et al., 2015) and the broader thermal range of pyroxene metastability (Bina, 2013) as shown in Figure 8 due to slow rates of diffusion-controlled pyroxene-garnet transformation in the transition zone (Nishi et al., 2013; Van Mierlo et al., 2013). However, it has not been observed in the laboratory, and both the low abundance of pyroxene in the peridotitic bulk of the slab (<20% by volume (Ringwood, 1991) and the generally displacive nature of metastable pyroxene transitions (Dera et al., 2013) also disfavor this mechanism. The arguments that the earthquake has a high seismic efficiency and fast rupture velocity V r already advanced against thermal instability also disfavor a melt-related mechanism such as option 2) for the first two rupture segments. Although the third rupture segment ruptures somewhat more slowly, the radiation efficiency η R   of 0.42, calculated via the equation η R = 1 − 1 − V r / β 1 + V r / β (Venkataraman and Kanamori, 2004), where β is shear velocity at the focal depth (5.5 km/s, Kennett et al., 1995), is much larger than 0.02, typical radiation efficiency of the thermal-induced earthquakes (Prieto et al., 2013). This precludes option 2) as the mechanism of the third rupture segment. Option 3) remains speculative because the 2013 Okhotsk and three prior transition-zone earthquakes reveal the same down-dip compressional stress regime in their focal mechanisms as that governing the intermediate-depth earthquakes (Figure 9). There is also no evidence of slab contortions in the seismicity, which displays a steadily dipping Wadati-Benioff zone (Figure 9), and seismic tomography suggests that the subducting Okhotsk slab proceeds into the lower mantle with roughly the same dip orientation (Fukao and Obayashi, 2013).

We consider the simplest plausible explanation to be that dehydration embrittlement (Green and Houston, 1995; Kirby, 1995; Houston, 2015) remains operable under conditions of deep-focus seismicity (Jung et al., 2004; Zhang et al., 2004). Jung et al. (2004) show that faulting due to dehydration can occur even at pressures so high that the volume change of reaction becomes negative and this mechanism can attain a near-critical state susceptible to stress change. This indicates that dehydration embrittlement does not require the evolution of a free fluid with a positive volume change to allow slip to occur. Ferrand et al. (2017) find that faults are experimentally produced by stress transfer driven by volume change of dehydration even under a partially hydrated condition (with only 5% volume of antigorite). Moreover, Green et al. (2015) and Green (2017) argue that it is nanometric crystallization of product phases, whether from hydration-dehydration reactions or from metastable transformations, which induce the formation of a fault zone. Deep dehydration has been invoked to interpret the rupture mechanism of the November 24, 2015 Mw 7.6 Peru deep-focus earthquake doublet (Zahradník et al., 2017).

Rüpke et al. (2006) demonstrate that at least 10% of water contained within the oceanic lithosphere upon subduction could survive passage through the mantle wedge en route to the mantle transition zone and lower mantle. While common hydrous minerals such as serpentine and chlorite would break down in the uppermost mantle, dense hydrous magnesium silicates can carry the water deeper and are stable in cool slabs down to at least as far as the top of the lower mantle (Komabayashi, 2006). Once these break down in/near the transition zone, the nominally anhydrous olivine polymorphs in the ambient mantle can absorb this water without melting (Smyth, 1994; Smyth and Jacobsen, 2006, studies in: Jacobsen and Van der Lee, 2006). Moreover, dehydration reactions in hydrated subducting slabs closely correlate with the depth-distribution of subduction-zone earthquakes, indicating dehydration embrittlement as a universal mechanism for intermediate and deep earthquakes (Omori et al., 2004; Barcheck et al., 2012; Shirey et al., 2019). In addition, there is evidence from deep diamond inclusions (Pearson et al., 2014) and indirect evidence from magneto-telluric studies, seismic tomography, and studies of seismic discontinuities (studies in Jacobsen and Van der Lee, 2006) that water indeed resides and cycles through the transition zone. Omori and Komabayashi (2007) propose one of several pathways by which subducting lithospheric slabs bring some of all subducted water to depth beyond the mantle wedge in dense hydrous magnesium silicates, which transform and dehydrate in the deep transition zone to help form less hydrous ringwoodite. Van der Lee et al. (2008) and Wang et al. (2018) show that hydrous ringwoodite is consistent with lowered rigidity imaged by seismic tomography in some places. The amounts by which the rigidity and density are lowered in ringwoodite are consistent with the volume of water supplied by a slab traversing the transition zone. Van der Lee et al. (2008) propose that a less-dense, hydrated transition zone provides a mechanism for water to return to the surface and that this deep water cycle may be responsible for many present and past plate-tectonic events. Shirey et al. (2019) invoke a similar deep water cycle to explain the deep diamond inclusions of Pearson et al. (2014). In other words, a deep upper-mantle water cycle is more consistent with diverse geophysical observations than a dry deep upper mantle.

Furthermore, Silver et al. (1995) suggested that the subhorizontal fault plane of the deep 1994 Mw 8.3 Bolivia earthquake could represent the reactivation of a fault that formed at the oceanic outer rise through bending of the subducting lithosphere, as is also argued for intermediate-depth seismicity (Jiao et al., 2000; Chen et al., 2004; Kiser et al., 2011). Warren et al. (2015) associated many such deep earthquakes with similar subhorizontal fault planes or their conjugate fault planes. Faulting at the oceanic outer rise has also been identified as an effective mechanism to hydrate the deep crust and mantle of the lithosphere before subduction (Rüpke et al., 2006). Moreover, due to the hydrated outer-rise normal fault, the amount of the subducted water is assessed to be one order of magnitude higher than previous estimation (Garth and Rietbrock, 2014; Cai et al., 2018), making these subhorizontal fault planes and their conjugates prime candidates for dehydration reactions and related embrittlement-induced and earthquake rupture.

If indeed dehydration played a role in the deep Okhotsk earthquake, the entirety of the subducting Pacific lithosphere cannot be completely dry even at these great depths in the mantle transition zone.

The 2013 Mw 8.6 Sea of Okhotsk deep-focus earthquake is found to rupture at different speeds to the north and the south, respectively, of the epicenter. This rupture model is constrained by a multi-array and multi-frequency combined back-projection method as well as a global P-wave train inversion technique for discrete subevents. The earthquake asymmetrically ruptures bilaterally a subhorizontal curved fault (150 km long, 70 km wide) in three stages. In the first 13 s, the rupture extends northwestwards on the downdip portion of the fault for approximately 55 km at a velocity of approximately 4.3 km/s. Afterward, a subevent at approximately 44 km southwest of the epicenter may be triggered by the passage of S-waves from previous faulting in the northwest and rupture continues eastward at a speed of 4.6 km/s. At the last stage, the rupture propagates up dip in a southeastward direction at a velocity of approximately 2.7 km/s until it ends at approximately 85 km SSE of the epicenter. The high rupture speeds indicate high seismic radiation efficiencies during the earthquake. Moreover, sizes of the subevents in the three rupture segments are beyond those of the possible MOW (Figure 8), which are not distorted as focal mechanisms of the nearby deep-focus earthquakes are subject to the similar down-dip compressional stress regime and seismicity is distributed regularly in the subducting slab. Thus, these findings together with diverse evidence of water subducting as deep as the transition zone and below suggest dehydration as the mechanism of the deep-focus earthquake.

• The 2013 M_w 8.3 Sea of Okhotsk deep-focus Earthquake bilaterally ruptured over 160 km in 30 s in three stages.