Deep Dehydration as a Plausible Mechanism of the 2013 Mw 8.3 Sea of Okhotsk Deep-Focus Earthquake

The rupture mechanisms of deep-focus (>300 km) earthquakes in subducting slabs of oceanic lithosphere are not well understood and different from brittle failure associated with shallow (<70 km) earthquakes. Here, we argue that dehydration embrittlement, often invoked as a mechanism for intermediate-depth earthquakes, is a plausible alternative model for this deep earthquake. Our argument is based upon the orientation and size of the plane that ruptured during the deep, 2013 Mw 8.3 Sea of Okhotsk earthquake, its rupture velocity and radiation efficiency, as well as diverse evidence of water subducting as deep as the transition zone and below. The rupture process of this earthquake has been inferred from back-projecting dual-band seismograms recorded at hundreds of seismic stations in North America and Europe, as well as by fitting P-wave trains recorded at dozens of globally distributed stations. If our inferences are correct, the entirety of the subducting Pacific lithosphere cannot be completely dry at deep, transition-zone depths, and other deep-focus earthquakes may also be associated with deep dehydration reactions.


INTRODUCTION
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 multiarray 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 ruptureassessment 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.

DATA AND METHODOLOGY
Multi-Array, Multi-Frequency-Band Back-Projection Analysis P-wave trains from two seismic networks TA and EU were back-projected to the source region (Zhang et al., 2016; in two adjacent frequency bands (0.1-0.5 Hz and 0.5-2.0 Hz, Figure 1). We developed the final back-projection rupture model in the following three steps.
The first step was to separately generate single-array backprojection imaging for each frequency band and each array. Because signal coherence is a critical factor in successful back projection, we only retained traces having a cross-correlation coefficient greater than 0.6 with the trace of a reference station near the array's center. This culling reduced the number of used traces to 74 high-frequency and 377 low-frequency traces for the TA array ( Figures 1B,E). Likewise, for the EU array, we retained 57 high-frequency and 69 low-frequency traces ( Figures 1C,F). The low-frequency data at the TA have better azimuthal coverage than the high-frequency data, whereas the azimuthal coverage for the EU is independent of the frequency band ( Figures 1A,D).
The global Centroid Moment Tensor (GCMT) implies two possible fault planes: one (strike 189 o , dip 11 o , and rake −93 o ) is subhorizontal, and the other one (strike 12 o , dip 79 o , and rake −89 o ) is subvertical ( Figures 1A,D). As elucidated by Meng et al. (2014), a subhorizontal plane agrees better with both the P and pP-waves back-projection imaging as well as our backprojection models. Moreover, the subhorizontal plane as the causative fault is evidently supported by a synthetic test done with two seismic point sources separating 20 km vertically as detailed in A Synthetic Test for the Causative Fault's Orientation and BP Images' Resolution. We approximated this subhorizontal fault plane by a horizontal plane centered around the USGS hypocenter (54.87 o N, 153.28 o E, and 609 km). The plane is larger than the area that contains the aftershocks (Supplementary Figure 1) and is gridded into 101 × 101 0.05 o × 0.05 o fault patches. A grid search was conducted to identify slip on those fault segments that correspond to stacking slownesses that maximize stacked P wave energy in 6-s-long time windows for the first 50 s of P arrivals (Figure 1).
In a second step, and for both frequency bands, the two sets of P-wave trains from the TA and EU arrays were jointly backprojected together. To equalize the contribution of the two arrays, stacked P amplitude distributions were normalized with their respective maxima. The multi-array back-projection imaging was Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 521220 then done by summing up the normalized single-array distributions and identifying the subevents for both frequency bands. The third step back-projected the high-frequency and lowfrequency data sets simultaneously. As before, the stacked P amplitude distributions for the two frequency bands were normalized with their respective maxima before stacking. Thus, the four sets of P-wave trains, from the TA and EU arrays and both frequency bands, were back-projected together to image slip evolution on the fault plane and in time ( Figure 3). We eliminated the artifacts from the subevents at the same locations by retaining the subevents with the maximum energy. The rupture time of a subevent was obtained by averaging those derived from both arrays. Thus, the remaining subevents comprise the final back-projection rupture model as shown in Figure 4. Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 521220 4

Global, Low-Frequency Waveform Inversion
Low-passed P-wave trains (2-100 mHz) from the Global Seismic Network (GSN) were inverted for seismic moments of the Okhotsk earthquake's subevents, using the method of Kikuchi and Kanamori (1991). The P waves are from 69 GSN stations with epicentral distances in the teleseismic window for a deep source from 20 o to 80 o (Figure 2A), avoiding involvement of P wave triplications and P waves in shadow zone. P-wave peak values in traces along the strike do not vary significantly with time ( Figure 2B), implying that all subevents may have similar nodal planes and thus focal mechanisms. In addition, the shortest duration of the P waveforms in the azimuthal range of 120 o -180 o represents rupture directivity associated with the southward rupture of the earthquake.
The subhorizontal fault plane as the causative fault, indicated by the BP synthetic test detailed in A Synthetic Test for the Causative Fault's Orientation and BP Images' Resolution, was gridded into 19 × 13 patches, each 10 km long along the strike and 5 km wide, respectively. Based on the subevents imaged in the final back-projection model, we chose to subdivide the rupture and slip evolution into a number of discrete subevents and found that a source model with 12 subevents (Supplementary Figure 1) produces well-fitting waveforms except for those with azimuths close to the strikes of nodal planes, such as station COLA (Supplementary Figure 2), and captures the complexity of the rupture pattern. Increasing the number of subevents improves the waveform fits but yields unrealistically high seismic moments (Supplementary Figure 3). To be consistent with the final backprojection rupture model, we removed subevents 40 s after the origin time, leading to a well correspondence of the remaining subevents to those in the multi-array, multi-frequency-band back-projection rupture model, except for three subevents northwest of the epicenter and the southernmost subevent ( Figure 4 and Supplementary Figure 1). Because the contribution of these four subevents to the waveform fits is very small (Supplementary Figure 2, 4), we removed them from the final multi-subevent model ( Figure 5A). The resulting model has a seismic moment of 2.3 × 10 21 Nm (equivalent to a moment magnitude of M w 8.2), providing a lower-bound for the seismic moment of the earthquake.
This multi-subevent model estimates moment release of the earthquake's subevents. Rupture areas of the subevents in the slip model can be inferred from their seismic moments. The seismic moment is defined as M 0 μuA (Aki, 1966), where μ is the shear modulus at the focal depth, u is the slip, and A is the rupture area. Assuming subevents rupture circularly gives A πR 2 , where R is the radius of the ruptured fault patch. For the earthquake, the average stress drop Δσ is estimated to be 12-15 MPa (Ye et al., 2013), which is compatible with the stress drop in a range of 5-50 MPa for 14 large deep earthquakes (Tibi et al., 2003). For a circular rupture, the rupture radius can then be estimated from the seismic moment as follows (Eshelby, 1957): The estimated rupture radii of the subevents range from 15.0 to 27.1 km ( Figure 5A). Estimating the rupture radii is important for evaluating the fault dimensions and inferring the source mechanism of the earthquake.

A SYNTHETIC TEST FOR THE CAUSATIVE FAULT'S ORIENTATION AND BP IMAGES' RESOLUTION
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  Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 521220 6 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 FIGURE 7 | Rupture imaging of the two seismic sources by back-projecting synthesized P-waves from the TA (left panels) and the EU (middle panels) in the lowfrequency band and the combined imaging (right panels). The color bar denotes the normalized power of potential subevents. Red stars indicate the epicenters. The contour shows the 80% of the maximum normalized power.
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 521220 7 from the single-array BP. This is the reason why we take multiarray, multi-band, BP imaging in this study.

RESULTS
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, multifrequency 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 highfrequency 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 lowfrequency 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 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 13secend 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).

COMPARISON WITH PREVIOUS STUDIES
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 singlearray 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 eventarray azimuth (Zhang et al., 2016;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.

DISCUSSION
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 intermediatedepth 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+Vr/β (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 deepfocus 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 deepfocus 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 andJacobsen, 2006, studies in: Jacobsen andVan 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 platetectonic 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 intermediatedepth 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.

CONCLUSION
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 downdip 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.

KEY POINTS
• The 2013 M w 8.3 Sea of Okhotsk deep-focus Earthquake bilaterally ruptured over 160 km in 30 s in three stages. • The deep earthquake ruptured fast and its rupture area extends beyond the inferred metastable olivine wedge. • The deep earthquake's rupture mechanism is consistent with it being the result of dehydration of deeply subducted oceanic lithosphere.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.