## Delineating Deep Faults

Most large, damaging earthquakes initiate in Earth's crust where friction and brittle fracture control the release of energy. Strong earthquakes can occur in the mantle too, but their rupture dynamics are difficult to determine because higher temperatures and pressures play a more important role. **Ye et al.** (p. 1380) analyzed seismic

*P*waves generated by the 2013

*M*

_{w}8.3 Sea of Okhotsk earthquake—the largest deep earthquake recorded to date—and its associated aftershocks. The earthquake ruptured along a fault over 180-kilometer-long and structural heterogeneity resulted in a massive release of stress from the subducting slab. In a set of complementary laboratory deformation experiments,

**Schubnel**(p. 1377) simulated the nucleation of acoustic emission events that resemble deep earthquakes. These events are caused by an instantaneous phase transition from olivine to spinel, which would occur at the same depth and result in large stress releases observed for other deep earthquakes.

*et al.*## Abstract

Earth’s deepest earthquakes occur in subducting oceanic lithosphere, where temperatures are lower than in ambient mantle. On 24 May 2013, a magnitude 8.3 earthquake ruptured a 180-kilometer-long fault within the subducting Pacific plate about 609 kilometers below the Sea of Okhotsk. Global seismic *P* wave recordings indicate a radiated seismic energy of ~1.5 × 10^{17} joules. A rupture velocity of ~4.0 to 4.5 kilometers/second is determined by back-projection of short-period *P* waves, and the fault width is constrained to give static stress drop estimates (~12 to 15 megapascals) compatible with theoretical radiation efficiency for crack models. A nearby aftershock had a stress drop one to two orders of magnitude higher, indicating large stress heterogeneity in the deep slab, and plausibly within the rupture process of the great event.

The occurrence of earthquakes in the depth range from 400 to 720 km (the mantle transition zone) has long been enigmatic, given the immense pressure exerted by the overlying rock mass on any fault surface. Seismic radiation from deep earthquakes indicates that they likely involve shear faulting basically indistinguishable from shallow earthquakes despite the extreme pressure conditions. Deep earthquakes only initiate in relatively low-temperature regions of subducted oceanic lithosphere. Very high deviatoric stresses may be present in the core of the subducted slab, and some mechanism must exist to offset the inhibiting effects of pressure to allow shear faulting to initiate (*1*). For the depth range from 50 to 400 km, it is generally believed that release of water by mineral dehydration reactions or production of other fluid phases reduces the effective normal stress on surfaces and enables fluid-assisted frictional sliding. It is not clear whether such mechanisms can account for transition-zone earthquakes, the largest of which tend to occur below a 600-km depth. Much research has focused on processes such as abrupt phase transitions (*2*) that may be able to nucleate rupture under tremendous confining stress. Once deep fault slip initiates and becomes substantial, frictional heating can lead to melting of the fault surface, abetting runaway rupture expansion for large deep earthquakes (*3*).

On 24 May 2013 the largest deep earthquake yet recorded occurred near a depth of 609 km [05:44:49 UTC, 54.874°N, 153.281°E (*4*)] in the Pacific plate subducting along the Kuril-Kamchatka subduction zone (Fig. 1). The event locates under the Sea of Okhotsk. Globally recorded long-period seismic waves indicate that the overall earthquake process appears to involve shear faulting with a seismic moment of ~4.1 × 10^{21} N·m [moment magnitude (*M*_{w}) = 8.3] (*4*, *5*). The event is slightly larger than the 637-km-deep Bolivia earthquake of 9 June 1994 that had a seismic moment of ~3 × 10^{21} N·m (*M*_{w} = 8.3) (*6*, *7*).

Both great events have similar faulting geometries with very shallow-dipping normal fault mechanisms and only minor deviations from shear double-couple solutions. The 1994 Bolivia earthquake was interpreted as having a relatively slow rupture velocity, *V*_{r} ~ 1 to 2 km/s, with an ~40-s rupture duration and a spatially compact rupture zone with a scale of about 40 km by 60 km (*6*–*11*), leading to large stress drop estimates of around 110 to 150 MPa (*7*, *11*).

In the first 4 days after the 2013 earthquake, nine aftershocks were detected, eight having small magnitudes of 4.1 to 4.4 at depths from 487 to 627 km, and an *M*_{w} = 6.7 event struck 9 hours after the mainshock [14:56:31, 52.222°N, 151.515°E, 623 km deep (*4*)] about 200 km to the southwest (Fig. 1). Six nearby aftershocks with depths of around 600 km define a north-south trend about 90 km long, preferentially extending southward from the mainshock hypocenter. The trend is generally compatible with rupture along either of the two nodal planes of the mainshock focal mechanism but slightly favors the shallow plane. Aftershock occurrence for large deep earthquakes is highly variable (*12*), and the large depth range for the 2013 aftershocks suggests that some are triggered away from the mainshock rupture zone. The 2013 event was preceded by nearby large earthquakes in 2008 (*M*_{w} = 7.3 and 7.7, Fig. 1), the larger of which was ~100 km along strike to the south.

Two of the most fundamental seismological properties of a large earthquake are the source spectrum and the radiated seismic energy. We analyzed extensive global seismic recordings of *P* waves for the 2013 event to estimate both. Source spectrum estimates obtained from two distinct approaches provide estimates of the radiated seismic energy (Fig. 2). The radiated energy estimates depend on attenuation corrections. The attenuation corrections are parameterized by *t**, the ratio of total *P* wave travel time to average attenuation quality factor, *Q*(*f*), along each path as a function of frequency, *f*. The precise *Q*(*f*) on each path is not known in detail and is expected to vary strongly because of upper mantle heterogeneity in attenuation structure beneath the seismic stations. For deep focus earthquakes, the *t** values are expected to be on average ~0.5 s at 0.1 Hz for teleseismic *P* waves (half of the *t** value for a shallow source) and about 0.25 to 0.5 s at 1.0 Hz, with *t** likely decreasing as frequency increases above 1.0 Hz. Lacking knowledge of specific path or even best average attenuation parameters, we show spectral estimates for a range of constant *t** values from 0.1 to 0.5 s, with a value of 0.3 s deemed to be a reasonable value. The uncertainty in *t** affects the high-frequency spectral levels, which are very important for radiated seismic energy estimates. We averaged the energy values obtained from 102 stations by integrating the energy spectrum from the *P*-wave ground-motion velocities (*13*–*16*) after correcting for faulting geometry and propagation effects. Use of *t** = 0.3 s for each station and frequencies up to 3 Hz gives radiated energy of *E*_{R} = 1.5 × 10^{17} J, with a range of reasonable estimates being given by results for *t** = 0.2 s (*E*_{R} = 1.0 × 10^{17} J) to 0.4 s (*E*_{R} = 2. 8 × 10^{17} J). For *t** = 0.3 s, the moment-scaled energy is *E*_{R}/*M*_{0} = 3.7 × 10^{−5}.

To confirm the source spectrum estimate, we used *P*-wave observations of the nearby *M*_{w} 6.7 aftershock at the same stations as for the mainshock to explicitly cancel out the unknown path effects. The large aftershock is remarkably short duration, with impulsive *P* wave motions that have average pulse widths of about 1.8 s. The average source spectrum for the aftershock found assuming *t** = 0.3 s for 22 stations has a very flat spectrum up to ~0.5 Hz (Fig. 2), indicating that this event can serve as an impulse response, or empirical Green’s function (EGF), up to near that frequency. For *t** = 0.3 s, *E*_{R} is 2.36 × 10^{15} J for the EGF from log averaging of the 22 individual station energies, and *E*_{R}/*M*_{0} is 2.8 × 10^{−4} using our finite-fault inversion estimates of the seismic moment *M*_{0} = 8.4 × 10^{18} N·m.

We computed mainshock/EGF spectral ratios for 284 broadband *P*-wave observations (fig. S1), correcting for differences in radiation pattern, geometric spreading, and multiplying by the EGF *M*_{0}. The ratios are in close agreement with the averages of mainshock *P* spectra over the corresponding passband (Fig. 2). This ensures that uncertainties in *t** do not bias the average spectrum estimate in this passband.

Extrapolations of the spectral ratios from 0.25 to 3.0 Hz are made for various assumed mainshock spectral decay slopes with frequency exponents varying from –1 to –2. For reference, a shallow interplate earthquake source spectrum for a moment equal to that of the mainshock has a decay exponent of –2, a stress parameter of 3 MPa, and source velocity of 3.75 km/s. The deep earthquake spectral amplitudes are enriched in high frequency relative to the reference model, in part because of higher source velocity and in part because of higher stress drop. The precise spectral decay slope expected near 1 Hz is not known, because it depends on the detailed space-time history of slip on the fault, but values around –1.5 to –2 are consistent with assuming *t** values around 0.3 s. We conclude that *E*_{R} ~ 1.5 × 10^{17} J, with about a factor of 2 uncertainty. This is about three times as large as for the 1994 Bolivia event [*E*_{R} ~ 5.2 × 10^{16} J (*3*, *17*)].

The spatial extent of the 2013 Sea of Okhotsk deep earthquake faulting is critical for estimating additional properties of the source, such as slip pattern and static stress drop. Back-projection of teleseismic short-period *P* waves was used to estimate *V*_{r} and the source rupture dimensions (Fig. 3). The data are from large continental seismic networks in Europe and North America that recorded coherent broadband waveforms (fig. S4). For both array geometries, the back-projections indicate asymmetric bilateral extent of short-period radiation extending 50 to 60 km to the north of the hypocenter and about 120 km to the south, along the trend of the deep aftershocks, with a source time duration of about 30 s. Animations of the back-projections show the space time evolution of the short-period radiation (movie S1).

If we adopt the *V*_{r} estimate of 4.0 km/s obtained from the back-projections as a constraint on the finite-fault inversions, we find rupture models with average slip of 1.9 to 2.3 m and an average static stress drop of 4 to 5 MPa for rupture areas that have a radius of about 74 to 82 km (figs. S2 and S3) for the two fault plane choices. For these estimates, only subfaults with moment at least 15% of the peak subfault moment were retained to diminish sensitivity to poorly resolved low-slip areas of the model (*18*).

A problem with these solutions is that they can give large (>1) estimates of calculated radiation efficiency, η_{R}, which is the ratio of *E*_{R} to the available potential energy, Δ*W*_{0}:η_{R} = *E*_{R}/Δ*W*_{0} = 2μ*E*_{R}/(Δσ_{s}*M*_{0})(1)where Δσ_{s} is the static stress drop and μ is the rigidity.

Radiation efficiency has been calculated as a function of *V*_{r} for mode II, mode III, and energy-based (mode E) crack models (Fig. 4) (*19*–*23*). For higher *V*_{r}, there is less energy dissipation near the crack tip, so the radiation efficiency approaches 1 as *V*_{r} approaches the limiting speed (the Rayleigh velocity for mode II and the shear velocity for mode III and mode E). For the 2013 Sea of Okhotsk event, *V*_{r} < 2.5 km/s is required for the circularly expanding rupture models to have large enough calculated stress drop to lower the seismic efficiency to intersect the predictions of crack theory (Fig. 4). Such a low *V*_{r} cannot account for the faulting dimensions indicated by the back-projections for the 2013 Sea of Okhotsk event. In order to obtain radiation efficiency consistent with the crack models, the width of the ruptures for high *V*_{r} must be constrained, essentially increasing the static stress drop by imposing rectangular rather than circular fault expansion.

For *V*_{r} = 4.0 km/s, Δσ_{s} = 15 MPa is needed to give η_{R} = 0.6 for a mode III rupture (Fig. 4). For a 180-km rupture length, imposing a fault width of 60 km yields an effective rupture area that gives the required stress drop. For *V*_{r} = 5 km/s, the fault width is increased to 68 km and gives Δσ_{s} = 12 MPa and η_{R} = 0.75 (Fig. 4). These models are now physically acceptable, and the slip distribution for the shallow-dipping plane for *V*_{r} = 4 km/s has average slip of 4.4 m (Fig. 1 and fig. S9a). Good fits to observed *P* waveforms are obtained (fig. S10). The slip distribution has asymmetric bilateral extent of 60 km north northeast and 120 km south southwest. Large slip is concentrated between 30 and 90 km south of the hypocenter, with the area of significant slip being 9675 km^{2}. We have some preference for the shallow-dipping plane, but very similar results are found for the steeply dipping plane (fig. S9b); the narrow rectangular faults give similar waveforms at most stations for the along-strike rupture.

Although there are limitations in the precision of the back-projection constraints and the theoretical crack-model efficiency is calculated for very simple models, the basic model with *V*_{r} ~ 4.0 km/s and Δσ_{s} ~15 MPa appears to be a valid representation of the overall source rupture. The *V*_{r} and rupture area are both about a factor of 4 larger than for the 1994 Bolivia event, and Δσ_{s} is about an order of magnitude lower. The Δσ_{s} ~15 MPa estimate is comparable to values for shallow intraplate earthquakes (*13*, *24*, *25*), and the fault dimensions are similar to those for the shallow trench-slope intraplate normal faulting event in the Kuril Islands of 13 January 2007 (*M*_{w} = 8.1) (*26*). The fault geometry is compatible with rerupture of such an outer rise normal fault surface within the plate, rotated by the dip of the deep slab. However, the *M*_{w} = 6.7 aftershock has an unusually short rupture duration and finite-fault inversions for variable assumed rupture velocities for that event give Δσ_{s} estimates in the range 157 to 5856 MPa (fig. S12). Independent constraints on the aftershock fault area or rupture velocity are not available, but there is no question that it has a localized stress drop greatly exceeding the average stress drop for the 8.3 mainshock and substantially lower radiation efficiency suggestive of a more dissipative source process. It is plausible that within the mainshock rupture zone there were corresponding very high stress drop slip patches that cannot be resolved. The envelope of teleseismic *P*-wave ground accelerations for the mainshock follows the source-time function shape (fig. S13), so a very heterogeneous stress distribution on the fault does appear likely (*27*–*29*), and the average parameters do not represent the total degree of slip and stress heterogeneity.

The 2013 mainshock rupture extends along the slab strike, with slip likely confined to the low temperature core of the slab. The subducted plate is older and colder than the slab where the 1994 Bolivia earthquake occurred, and this difference in thermal state may have fundamentally affected how rupture expanded for the two events (*30*). The Sea of Okhotsk event is similar to a shallow intraplate earthquake, with a large aspect ratio fault area defined by the brittle core of the slab. For the Bolivia event, the brittle core volume is smaller, and surrounding ductile or plastic material with a finite strength may dominate. Faulting for the Bolivia event involved a rupture with a very dissipative source process that deposited a large amount of energy into the rupture zone, likely leading to melting. This behavior may be akin to that of a shear band. The Sea of Okhotsk mainshock rupture appears to have been less dissipative, and little or no melting may have occurred, although seismology cannot directly constrain the amount of melting.

The stress drops found for the 2013 mainshock and large aftershock suggest preexisting zones with strong and weak regions, likely inherited from shallow faulting of the slab. The warm plate in Bolivia may have more strong, less-brittle patches than weak-brittle patches, whereas the cold plate in the Sea of Okhotsk has more weak-brittle patches than strong, less-brittle patches. Strong patches may be distributed only sparsely in the Sea of Okhotsk slab, with one rupturing in the aftershock, and may act to stop rupture propagation on weak patches. The difference in the distribution of strong, less-brittle and weak-brittle patches caused by the difference in the thermal state is likely responsible for the drastically different source characteristics of the Bolivian and the Sea of Okhotsk events.

## Supplementary Materials

www.sciencemag.org/content/341/6152/1380/suppl/DC1

Materials and Methods

Figs. S1 to S13

References

Movie S1

## References and Notes

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
**Acknowledgments:**We thank E. Brodsky for helpful discussions and three anonymous reviewers for their thoughtful comments. The Incorporated Research Institutions for Seismology (IRIS) data management center provided the seismic recordings. This work was supported in part by NSF under grant EAR-1245717 (T.L.). All data used are available from the IRIS data center at www.iris.edu.