Energy Release of the 2013 Mw 8.3 Sea of Okhotsk Earthquake and Deep Slab Stress Heterogeneity

See allHide authors and affiliations

Science  20 Sep 2013:
Vol. 341, Issue 6152, pp. 1380-1384
DOI: 10.1126/science.1242032

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 Mw 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 et al. (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.


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 × 1017 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 × 1021 N·m [moment magnitude (Mw) = 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 × 1021 N·m (Mw = 8.3) (6, 7).

Fig. 1 Tectonic setting of the 2013 Mw 8.3 deep Sea of Okhotsk slab earthquake.

(Inset) The plate configuration, with the Pacific plate underthrusting the North American/Sea of Okhotsk plate along the Kuril-Kamchatka subduction zone at a convergence velocity of ~8.0 cm/year. Dashed lines are depth contours for the subducted oceanic slab beneath the Sea of Okhotsk. The main map shows best-double couple faulting mechanisms from global centroid-moment tensor inversions for recent large earthquakes in the deep slab, with blue indicating events below 600 km depth and cyan indicating events around 500 km deep. Focal mechanisms are at the event epicenter unless offset with a tie line. Small circles are locations of aftershocks of the 24 May 2013 event with magnitudes from 4.1 to 4.4. The contoured plot indicates the slip distribution of the preferred rupture model for the mainshock, with the large red star being the hypocenter at a depth of 608.9 km. The arrows indicate the magnitude and the direction of slip of the upper side of the fault, with the fault dipping 10° toward the west. The peak slip is 9.9 m, and the colors indicate about 2.2-m slip contour increments.

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, Vr ~ 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 (611), 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 Mw = 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 (Mw = 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 (1316) 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 ER = 1.5 × 1017 J, with a range of reasonable estimates being given by results for t* = 0.2 s (ER = 1.0 × 1017 J) to 0.4 s (ER = 2. 8 × 1017 J). For t* = 0.3 s, the moment-scaled energy is ER/M0 = 3.7 × 10−5.

Fig. 2 Source spectra estimates for the 2013 Mw 8.3 Sea of Okhotsk mainshock and Mw 6.7 aftershock.

The mainshock spectrum is estimated by two methods. The blue curves are estimates based on the spectrum of the source time function from finite-fault inversion for frequencies less than 0.05 Hz and from averaging many teleseismic P wave spectra with propagation and radiation pattern corrections from 0.05 to 3.0 Hz. Results are shown for different attenuation parameters of t* = 0.1 to 0.5 s. The red curve is an estimate of the mainshock source spectrum from 284 spectral ratios of the mainshock and the aftershock (EGF) spectra at the same station for the frequency band 0.03 to 0.25 Hz. The extrapolated spectra to 3 Hz assume source spectrum decay exponents from –1.0 to –2.0. The green curve is the average source spectrum for the Mw 6.7 event based on the first method, using an assumed t* = 0.3 s. (Inset) The dependence of estimates of ER on the assumed value of t* for the mainshock signals, given by averaging energy estimates from individual path-corrected P wave spectra.

To confirm the source spectrum estimate, we used P-wave observations of the nearby Mw 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, ER is 2.36 × 1015 J for the EGF from log averaging of the 22 individual station energies, and ER/M0 is 2.8 × 10−4 using our finite-fault inversion estimates of the seismic moment M0 = 8.4 × 1018 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 M0. 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 ER ~ 1.5 × 1017 J, with about a factor of 2 uncertainty. This is about three times as large as for the 1994 Bolivia event [ER ~ 5.2 × 1016 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 Vr 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).

Fig. 3 Constraints on rupture velocity from P wave back-projection.

Teleseismic P waves in the frequency band from 0.5 to 2.0 Hz from large networks of stations in North America (NA) and Europe (EU) were used to image the space-time history of coherent high-frequency seismic radiation from the 2013 Mw 8.3 Sea of Okhotsk earthquake. The time-integrated power stacked on a grid around the source are shown here relative to the mainshock epicenter (white stars). The darker blue colors indicate coherent energy release with an asymmetric spread of source radiation in the north-south direction being resolved by the images from both networks. movie S1 shows the time-varying images throughout the rupture process.

If we adopt the Vr 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 ER to the available potential energy, ΔW0:ηR = ERW0 = 2μER/(ΔσsM0)(1)where Δσs is the static stress drop and μ is the rigidity.

Radiation efficiency has been calculated as a function of Vr for mode II, mode III, and energy-based (mode E) crack models (Fig. 4) (1923). For higher Vr, there is less energy dissipation near the crack tip, so the radiation efficiency approaches 1 as Vr 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, Vr < 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 Vr 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 Vr must be constrained, essentially increasing the static stress drop by imposing rectangular rather than circular fault expansion.

Fig. 4 Model constraints from consideration of radiation efficiency for crack models with varying rupture speed.

Reference curves for mode II and mode III cracks and an energy-based model (mode E) have ηR values that approach 1 as the rupture speeds approach their limiting velocities (~5.1 km/s for mode II, ~5.5 km/s for mode III and mode E at a depth of 610 km) (16, 23). The blue circles indicate calculated radiation efficiencies for rupture models (fig. S2) with varying Vr and fault dimensions scaling proportional to Vr. These models are only compatible with the crack theory for Vr ~1.5 to 2.0 km/s, but this is inconsistent with the rupture extent indicated by back-projection in Fig. 3, which favors Vr of 4.0 to 5.0 km/s. By constraining the width of the slip models, we find high-Vr models consistent with the theoretical radiation efficiency, with preferred models giving the red stars. The solution shown in Fig. 1 is the 15 MPa stress drop model for Vr = 4.0 km/s. The pink bars indicate variation in estimates for different thresholds (0.1 to 0.2, with stars for 0.15) used to remove poorly resolved low slip regions of the fault models.

For Vr = 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 Vr = 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 Vr = 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 km2. 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 Vr ~ 4.0 km/s and Δσs ~15 MPa appears to be a valid representation of the overall source rupture. The Vr 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 (Mw = 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 Mw = 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 (2729), 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

Materials and Methods

Figs. S1 to S13


Movie S1

References and Notes

  1. 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
View Abstract

Stay Connected to Science

Navigate This Article