## Abstract

Rupture properties, such as rupture direction, length, propagation speed and source duration, provide important insights into earthquake mechanisms. One approach to estimate these properties is to investigate the body-wave duration that depends upon the relative location of the station with respect to the rupture direction. Under the assumption that the propagation is unilateral, the duration can be expressed as a function of the dip and azimuth of the rupture. Examination of duration measurements with respect to both the take-off angle and the azimuth is crucial to obtain robust estimates of rupture parameters, especially for nearly vertical rupture propagation. Moreover, limited data coverage, such as using only teleseismic data, can bias the source duration estimate for dipping ruptures, and this bias can map into estimates of other source properties such as rupture extent and rupture speed. Based upon this framework, we introduce an inversion scheme that uses the duration measurements to obtain four parameters: the source duration, a measure of the rupture extent and speed, and dip and azimuth of the rupture propagation. The method is applied to two deep-focus events in the Sea of Okhotsk region, an Mw 7.7 event that occurred on 2012 August 14 and an Mw 8.3 event from 2013 May 24. The source durations are 26 ± 1 and 37 ± 1 s, and rupture speeds are 49 ± 4 per cent and 26 ± 3 per cent of shear wave speed for the Mw 7.7 and 8.3 events, respectively. The azimuths of the two ruptures are parallel to the trench, but are in opposite directions. The dips of the Mw 7.7 and 8.3 events are constrained to be 48° ± 8° downdip and 19° ± 8° updip, respectively. The fit to the data is significantly poorer for the Mw 8.3 event than the Mw 7.7 event, suggesting that the unilateral rupture may not be a good assumption. The analysis is expanded into a multi-episode model, and a secondary episode is determined for the Mw 8.3 event in the southeast direction. The two-episode model gives a better fit to the data than the unilateral model and is compatible with the back-projection analysis, demonstrating that the rupture propagation of the Mw 8.3 event is complex.

## INTRODUCTION

Characteristics of earthquakes such as rupture direction can be investigated using teleseismic recordings. For example, finite-fault modelling (e.g. Kikuchi & Kanamori 1991; Wald et al.1993; Lay et al.2010; Hayes 2011) and back-projection method (e.g. Ishii et al.2005; Xu et al.2009; Koper et al.2011; Meng et al.2011; Yagi et al.2012) use waveform information, and they show time evolution of slip or energy release locations. Another approach to resolve the earthquake source properties is to investigate the directivity effect, that is, differences in the observed P-wave duration (e.g. Haskell 1964; Beck et al.1995; Warren & Silver 2006; Frez et al.2010; Lengliné & Got 2011; Bezada & Humphreys 2012; Convertito et al.2012). Examination of the observed durations allows estimates of source duration, rupture speed and rupture direction to be made (Haskell 1964), but it can be challenging due to limited data coverage and/or quality (Warren & Shearer 2006; Caldeira et al.2010). To make the problem tractable, assumptions or a priori constraints are often imposed. One of the common assumptions is the horizontal rupture propagation, that is, the observed durations depend only on the station azimuth (e.g. Pino et al.1999; Pro et al.2007). Other studies assume that rupture occurs on pre-determined nodal planes (e.g. Caldeira et al.2010). Different assumptions can lead to different results for the same earthquake, and it becomes difficult to compare the results (Frohlich 2006).

In this paper, we demonstrate that the 3-D directivity effects can be modelled, and that the duration variation in the dip direction is as important as that in the azimuthal direction. Based on the framework, an inversion method to estimate rupture parameters is introduced, which requires no a priori assumptions about the geometry of the earthquake such as fault plane orientation or horizontal propagation. Four parameters, that is, source duration, a measure of rupture extent and speed, rupture dip and azimuth, are obtained for a unilateral rupture model, and their uncertainties are assessed. Moreover, a method to analyse complex, non-unilateral ruptures is presented. We apply the technique to two deep-focus events and compare the parameters against back-projection results.

## METHOD

We seek a 3-D expression relating observed duration to dip and azimuth of rupture propagation assuming unilateral rupture with a constant rupture speed. The formulation demonstrates that bias in a source duration estimate arises even when the azimuthal coverage is perfect but dip coverage is limited. Based upon this framework, a method combining grid search and iterative inversion is set up to obtain four rupture parameters: the source duration, a measure of rupture extent, dip and azimuth of the rupture propagation.

### Observed duration in 3-D

For a unilaterally rupturing earthquake, the observed duration, τi, at station i depends upon the relative location of the station with respect to the orientation of the rupture direction such that (Haskell 1964)

(1)
$$\tau _i = \tau -\frac{L}{V} \cos \theta _i,$$
where τ is the true source duration, that is, the time over which seismic energy is released, L is the rupture length, V is the seismic wave speed at the source and θi is the angle between the rupture and the take-off directions for the ith station (Fig. 1). The factor cos θi depends upon a combination of the dip and azimuth of the rupture and the ray path towards the station, that is,
(2)
$$\cos \theta _i = \sin \gamma \sin \gamma _i + \cos \gamma \cos \gamma _i \cos (\phi -\phi _i),$$
where (γ, ϕ) and (γi, ϕi) represent dip and azimuth of the rupture and the take-off directions, respectively. Here, dip is defined such that negative values correspond to updip rupture or waves taking off in the upper focal hemisphere. Combining the two expressions, the observed duration becomes
(3)
$$\tau _i = \tau -\frac{L}{V} \big [\sin \gamma \sin \gamma _i + \cos \gamma \cos \gamma _i \cos (\phi -\phi _i)\big ].$$
This explicit 3-D formulation allows one to better understand the directivity effect of a rupture with an arbitrary orientation. Furthermore, the source duration can be defined as the average over a unit sphere,
(4)
$$\bar{\tau } = \frac{1}{4\pi }\int _{-\frac{\pi }{2}}^{\frac{\pi }{2}} \int _0^{2\pi } \tau _i \cos \gamma _i \text{d} \phi _i \text{d} \gamma _i.$$
Substituting eq. (3) into the above equation recovers the source duration,
(5)
$$\bar{\tau } = \tau .$$

Figure 1.

Rupture direction, take-off vector and relative angles. Rupture direction vector (thick solid arrow) is shown on a shaded fault plane with its dip γ and azimuth ϕ. The take-off vector (dashed arrow) of seismic ray (thin black curve) to the ith station has dip γi and azimuth ϕi. The angle between the rupture direction and the take-off vector is θi.

Figure 1.

Rupture direction, take-off vector and relative angles. Rupture direction vector (thick solid arrow) is shown on a shaded fault plane with its dip γ and azimuth ϕ. The take-off vector (dashed arrow) of seismic ray (thin black curve) to the ith station has dip γi and azimuth ϕi. The angle between the rupture direction and the take-off vector is θi.

### Effect of rupture dip

A number of papers have used eq. (1) where θi is taken to be the azimuth of the station with respect to the earthquake (e.g. Pino et al.1999; Ammon et al.2005; Ni et al.2005; Pro et al.2007). This assumes that the rupture propagates on a horizontal plane, that is, dip of zero. However, ruptures are not always horizontal, and often, we do not have sufficient knowledge about the rupture direction or the fault plane, especially for intermediate or deep earthquakes (Persh & Houston 2004; Frohlich 2006). By examining the dependence of the observed duration on the dip angle, problems associated with horizontal rupture assumption or limited data coverage can be demonstrated.

To illustrate the issues, three cases, γ = 0°, 45° and 90°, are considered. In the horizontal rupture case (γ = 0°), the observed duration becomes

(6)
$$\tau _i = \tau -\frac{L}{V} \cos \gamma _i \cos (\phi - \phi _i).$$
The dependence of this equation on γi implies that even for a horizontal rupture propagation, the durations measured at stations with the same azimuth but at various distances, that is, various take-off angles, are different (Fig. 2a). Projection of the observed duration onto the focal sphere shows the same pattern in the upper and lower hemispheres (Fig. 3a), and the shortest and longest durations are observed for rays taking off in directions parallel and opposite to the rupture, respectively (Figs 2a,b and 3a). In contrast, when the rupture propagates in the vertical direction (γ = 90°), the observed duration depends only upon the dip of the take-off vector, that is,
(7)
$$\tau _i = \tau -\frac{L}{V} \sin \gamma _i,$$
such that there is no dependence on azimuth (Figs 2e,f and 3a). The lower-hemispheric projection of the duration shows that the shortest duration occurs at γi = 90° (Fig. 3a), implying that the use of teleseismic data (i.e. lower-hemispheric data) will result in a biased source duration estimate. The upper hemisphere has exactly the opposite pattern, and only the stations with waves taking off horizontally from the source will observe the true duration. In the final example where the rupture dip is 45°, the observed durations are
(8)
$$\tau _i = \tau -\frac{L}{V} \left[\frac{1}{\sqrt{2}}\sin \gamma _i+\frac{1}{\sqrt{2}}\cos \gamma _i \cos (\phi - \phi _i) \right].$$
In this case, the durations show features between those of the horizontal and vertical cases (Figs 2c,d and 3a). The shortest and longest durations occur at 45° from the horizontal plane.

Figure 2.

(a) Observed duration as a function of take-off dip (γi; Fig. 1) for horizontal rupture propagating northward (ϕ = 0°). The black and grey solid lines are for take-off azimuths of ϕi = 0° and ϕi = 180°, respectively, and the black and grey dashed lines are for ϕi = 45° and ϕi = 135°, respectively. The dotted line is when ϕi = 90°. (b) Observed duration as a function of take-off azimuth for horizontal rupture. The dotted line is for take-off dip of γi = 0°, the dashed line is for γi = ±45° and the solid line is for γi = ±90°. (c) Same as (a) except for a 45°-dipping rupture case. (d) Same as (b) except for a 45°-dipping rupture case. The black and grey solid lines are for γi = 90° and γi = −90°, respectively, and the black and grey dashed lines are for γi = 45° and γi = −45°, respectively. (e) Same as (a) except for a vertical rupture case. The dependence on the dip angle is the same for all azimuths. (f) Same as (d) except for a vertical rupture case.

Figure 2.

(a) Observed duration as a function of take-off dip (γi; Fig. 1) for horizontal rupture propagating northward (ϕ = 0°). The black and grey solid lines are for take-off azimuths of ϕi = 0° and ϕi = 180°, respectively, and the black and grey dashed lines are for ϕi = 45° and ϕi = 135°, respectively. The dotted line is when ϕi = 90°. (b) Observed duration as a function of take-off azimuth for horizontal rupture. The dotted line is for take-off dip of γi = 0°, the dashed line is for γi = ±45° and the solid line is for γi = ±90°. (c) Same as (a) except for a 45°-dipping rupture case. (d) Same as (b) except for a 45°-dipping rupture case. The black and grey solid lines are for γi = 90° and γi = −90°, respectively, and the black and grey dashed lines are for γi = 45° and γi = −45°, respectively. (e) Same as (a) except for a vertical rupture case. The dependence on the dip angle is the same for all azimuths. (f) Same as (d) except for a vertical rupture case.

Figure 3.

(a) Predicted duration of horizontal (left column), 45°-dipping (middle column) and vertical (right column) ruptures plotted on the upper (top row) and the lower (bottom row) focal hemispheres in the stereographic projection. Strike of the rupture is north (ϕ = 0°). Colours represent the observed duration. White solid lines are contour of the duration at 20 per cent interval, and the arrow head and tail notations in white indicate the rupture direction. (b) Bias resulting from lower-hemispheric data coverage. The difference between the estimated and source duration is shown on the y-axis.

Figure 3.

(a) Predicted duration of horizontal (left column), 45°-dipping (middle column) and vertical (right column) ruptures plotted on the upper (top row) and the lower (bottom row) focal hemispheres in the stereographic projection. Strike of the rupture is north (ϕ = 0°). Colours represent the observed duration. White solid lines are contour of the duration at 20 per cent interval, and the arrow head and tail notations in white indicate the rupture direction. (b) Bias resulting from lower-hemispheric data coverage. The difference between the estimated and source duration is shown on the y-axis.

These examples illustrate the problems of estimating the source duration using only the teleseismic observations, or, in the focal sphere projection, the lower-hemispheric data coverage. Based upon the definition of the average duration (eq. 4), the lower-hemispheric coverage gives

(9)
$$\bar{\tau }^l= \frac{1}{2\pi }\int _{0}^{\frac{\pi }{2}} \int _0^{2\pi } \tau _i \cos \gamma _i \text{d} \phi _i \text{d} \gamma _i = \tau -\frac{L}{2V} \sin \gamma ,$$
where $$\bar{\tau }^l$$ is the average using only the lower-hemispheric data. If the value is taken as the source duration estimate, the bias in the estimate depends sinusoidally on the dip of the rupture γ (Fig. 3b). It indicates that the source duration is overestimated for updip ruptures (γ < 0), and underestimated for downdip ruptures (γ > 0). As the rupture direction moves away from the horizontal plane, the deviation in the estimate from the true duration becomes larger. Furthermore, data coverage is uneven in real applications, and stations cover only a narrow range in dip (i.e. distance). In order to account for these problems, the dependence on dip should be included, for example, when teleseismic data are stacked to estimate the source duration or the source time function (e.g. Vidale & Houston 1993; Bos et al.1998; Houston et al.1998; Persh & Houston 2004).

### Inversion for the rupture properties

Based upon eq. (3), the duration measurements are modelled in terms of four parameters, source duration τ, ratio of rupture extent to the product of source duration and compressional wave speed k = LV, dip of the rupture γ and azimuth of the rupture ϕ. Note that the source duration, τ, is the sum of rupture duration, L/VR, and the rise time, where VR is the rupture speed. This implies that the value of the parameter k represents the ratio of rupture speed to compressional wave speed, VR/V, only if the rise time is zero, and if the rise time is finite, k gives an underestimated speed ratio. The rupture length L can be recovered from the value of k using wave speed at the source and we report both k and L estimates. Rewriting eq. (3) with the four parameters,

(10)
\begin{eqnarray} &&{\tau _i (\tau ,k,\gamma ,\phi )}\nonumber\\ &&{\quad= \tau \Big \lbrace 1-k\big [\sin \gamma \sin \gamma _i + \cos \gamma \cos \gamma _i \cos (\phi -\phi _i)\big ]\Big \rbrace .} \end{eqnarray}
For the inverse problem, we define the misfit function as
(11)
$$f(\tau ,k,\gamma ,\phi ) = \sum _{i=1}^N w_i \big [\tau _i (\tau ,k,\gamma ,\phi ) -\tau _i^{\text{obs}}\big ]^2,$$
where the summation is over N measurements, $$\tau _i^{\text{obs}}$$ is the observed duration for the ith station and wi is the weighting assigned to each observation based upon factors such as the quality of the measurement and station distribution. Since this is a nonlinear problem for the four parameters, a two-step procedure is introduced. The first step finds the solution for the parameters through grid search over rupture dip and azimuth angles. The second step uses the results from the first step as the starting model, and iteratively solves for perturbations in model parameters using a linearized form of eq. (10). The two steps are described in detail below.

The first step of obtaining the starting model parameters τ0, k0, γ0 and ϕ0 is achieved through a combination of analytical and grid search approaches. The initial source duration τ0 is calculated as the weighted average of the duration measurements,

(12)
$$\tau ^0 = \frac{\sum w_i \tau _i^{\text{obs}}}{\sum w_i},$$
where the summation over i from 1 to N is understood. Using this source duration, grid search over dip and azimuth is carried out, where the optimal k value is computed for each value of γ and ϕ using an expression obtained from minimizing the misfit function (eq. 11), that is,
(13)
\begin{eqnarray} k(\gamma , \phi ) &=& \frac{\sum w_i (\tau ^0 - \tau _i^{\text{obs}})\xi _i}{\tau ^0 \sum w_i \xi _i^2}, \quad \mathrm{where}\nonumber\\ \xi _i &=& \sin \gamma \sin \gamma _i + \cos \gamma \cos \gamma _i \cos (\phi -\phi _i). \end{eqnarray}
This parameter k is not constrained to be positive. As a result, the grid search returns two minima in misfit over the dip and azimuth space, corresponding to positive and negative k values. The two solutions are identical, that is, ruptures propagating in opposite directions with the velocities of the opposite sign. Therefore, the positive k value is taken to be k0.

Assuming that the grid search result is close to the optimal solution, eq. (10) is linearized to facilitate iterative inversion. The perturbation in observed duration δτi depends upon perturbations in source duration δτ, δk, dip δγ and azimuth δϕ as

(14)
\begin{eqnarray} \delta \tau _i &=& (1-k\xi _i) \delta \tau - \tau \xi _i \delta k \nonumber\\ &&-\,\tau k \left[ \cos \gamma \sin \gamma _i - \sin \gamma \cos \gamma _i \cos (\phi -\phi _i)\right] \delta \gamma \nonumber\\ &&+\, \tau k \cos \gamma \cos \gamma _i \sin (\phi -\phi _i) \delta \phi , \end{eqnarray}
after ignoring the second and higher order terms. Note that the factor multiplying δτ is dimensionless whereas the others are in units of time. In order to ensure consistency within the sensitivity matrix, the source duration term is non-dimensionalized by τ0, that is, τ0(1 − kξi) δ(τ/τ0). The non-dimensionalized expression of eq. (14) can be written in the matrix form,
(15)
$$\mathbf {\delta d}_{j-1} \, = \mathbf {G}_{j-1} \cdot \mathbf {\delta m} _j$$
where δdj − 1 is the vector containing the residual duration of N stations from the (j − 1)th iteration, [δτ1  δτ2 … δτN]T and δmj is the perturbation in the non-dimensionalized model parameters, [δτ/τ0  δk  δγ  δϕ]T, at the jth iteration. The elements of the Gj − 1 are calculated using the model parameters τ, k, γ and ϕ from the (j − 1)th iteration. In terms of the perturbation vector δd, and an N × N diagonal weighting matrix W with elements wi, the misfit function (eq. 11) becomes
(16)
$$f_j(\tau ,k,\gamma ,\phi ) = \mathbf {\delta d}_j^T \cdot \mathbf {W \cdot \delta d}_j.$$
Minimizing this weighted misfit function produces the model update,
(17)
$$\mathbf {\delta m}_j = (\mathbf {G}^T_{j-1} \cdot \mathbf {W \cdot G}_{j-1})^{-1} \cdot \mathbf {G}^T_{j-1} \cdot \mathbf {W \cdot \delta d}_{j-1},$$
with the updated model parameter vector,
(18)
$$\mathbf {m}_j = \mathbf {m}_{j-1} + \mathbf {\delta m}_j .$$
This inverse problem is over-determined and requires no regularization (except for some unusual circumstances such as when stations are all located on the nodal planes).

### Multi-episode directivity analysis for complex rupture

Thus far, the discussion assumed unilateral rupture, but in practice, there have been observations of complex rupture patterns (e.g. Butler et al.1979; Yagi & Kikuchi 2000; Zhang et al.2009; Hayes et al.2010). Let us consider a complex source where the rupture consists of multiple episodes occurring at arbitrary times and locations. These episodes can be large bursts of seismic energy or start/termination of slip, and the timing of the nth episode observed at the ith station, tni, is dependent not only on the episode's time tn, but also on its location, that is,

(19)
\begin{eqnarray} &&{t_{ni}(t_n,L_n,\gamma _n,\phi _n)} \nonumber\\ &&{\quad= t_n-\frac{L_n}{V}[\sin \gamma _n \sin \gamma _i + \cos \gamma _n\cos \gamma _i\cos (\phi _n-\phi _i)],} \end{eqnarray}
where Ln, γn and ϕn are the 3-D location of the nth episode, that is, distance, dip and azimuth, with respect to the hypocentre. Comparison of this expression with eq. (3) shows that the episode timing tn is the source duration τ in the unilateral case. Similarly, rupture properties estimated by the unilateral directivity analysis provides the timing and the location of the final episode for earthquakes with complicated rupture.

Eq. (19) also shows that the signal from the last episode is not always the last arrival at a given station. Depending upon the location of the station as well as location and timing of the episodes, waves generated by an earlier episode can be recorded after those from the final episode. For a given station i, the measured ‘end time’ is simply the latest of the episode arrivals, that is,

(20)
$$\tau _i = \max \left(t_{ni} \right).$$
The observed durations of an earthquake with complex rupture, therefore, do not necessarily follow a sinusoidal pattern. For example, when two episodes occur at the same time but at opposite sides of the hypocentre such as in bilateral rupture (Fig. 4a), long durations are observed at stations in the directions of the two episodes with respect to the hypocentre (Figs 4a–c). The shortest durations occur in a plane perpendicular to the line of rupture, and the observed durations do not vary smoothly but have cusps (Fig. 4b). If the two episodes are 90° from one another with respect to the hypocentre with the first episode located farther from the hypocentre yet happening earlier than the second episode (Fig. 4d), the longest duration is observed at stations that are farthest away from the two episodes (Figs 4e and f). In both examples, parts of the sinusoid associated with each episode are not observed (Figs 4b and e). The average duration (eq. 4) for a complex event, therefore, results in a value that is longer than the true source duration. This implies that applying the unilateral rupture analysis to a complex event leads to an overestimated source duration. Instead of the average, the source duration of a multi-episode event should be defined as the time between the occurrence of the last episode and the hypocentre. Also note that in the multi-episode rupture model, the shortest durations do not correspond to the rupture propagation directions or the episode locations, and it is more informative to search for the occurrences of longest durations (Fig. 4c).

Figure 4.

Examples of duration complexities for non-unilateral ruptures. (a) Depiction of a bilateral rupture on a horizontal plane. The hypothetical rupture starts at the white star (hypocentre), and propagates bilaterally to south (green star, episode 1) and north (blue star, episode 2). The two episodes are at an equal distance from the hypocentre and occur at the same time. (b) Episode timing (y-axis) as a function of take-off azimuth for the model depicted in (a). The dotted lines are for take-off dip of 0° with colours matching the corresponding episode, and the dashed lines are for take-off dip of ±45°. Thick red curves indicate the end times. (c) Event duration in the focal hemisphere stereographic projection. The pattern for the upper and lower hemispheres are identical. White solid curves are contours of the duration at 20 per cent interval, and the black dotted and dashed circles correspond to the dip angles shown in (b). (d) Same as (a), except that episode 1 (green star) occurs earlier than episode 2 (blue star), and at a distance farther from the hypocentre. (e) Same as (b) except that the lines are for the model shown in (d). (f) Same as (c) except that the plot is for the model shown in (d).

Figure 4.

Examples of duration complexities for non-unilateral ruptures. (a) Depiction of a bilateral rupture on a horizontal plane. The hypothetical rupture starts at the white star (hypocentre), and propagates bilaterally to south (green star, episode 1) and north (blue star, episode 2). The two episodes are at an equal distance from the hypocentre and occur at the same time. (b) Episode timing (y-axis) as a function of take-off azimuth for the model depicted in (a). The dotted lines are for take-off dip of 0° with colours matching the corresponding episode, and the dashed lines are for take-off dip of ±45°. Thick red curves indicate the end times. (c) Event duration in the focal hemisphere stereographic projection. The pattern for the upper and lower hemispheres are identical. White solid curves are contours of the duration at 20 per cent interval, and the black dotted and dashed circles correspond to the dip angles shown in (b). (d) Same as (a), except that episode 1 (green star) occurs earlier than episode 2 (blue star), and at a distance farther from the hypocentre. (e) Same as (b) except that the lines are for the model shown in (d). (f) Same as (c) except that the plot is for the model shown in (d).

The relationship between observed duration and episode parameters (eqs 19 and 20) can be used to obtain models of episode occurrences. The same misfit function given by eq. (11) can be used, but the optimization problem becomes nonlinear with non-differentiable points (Figs 4b and e). A global optimization scheme is needed, and we adopt the simulated annealing algorithm (Metropolis et al.1953; Kirkpatrick et al.1983; Ingber 1995). In order to ensure convergence toward the global minimum, one of the parameters, temperature which controls the probability of moving to a different set of model parameters, is reduced logarithmically, that is, the temperature at the jth iteration Tj is equal to T0/ln j, where T0 is the initial temperature (Geman & Geman 1984).

## DATA

The inversion scheme for the rupture directivity analysis is applied to two events in the Sea of Okhotsk region, an Mw 7.7 event that occurred on 2012 August 14 (Event 1) and an Mw 8.3 event from 2013 May 24 (Event 2). According to the USGS National Earthquake Information Center (NEIC), the hypocentres of Events 1 and 2 are located at 49.80°N, 145.06°E and 583-km depth, and at 54.87°N, 153.33°E and 601.7-km depth, respectively (Fig. 5). These deep-focus events are chosen, because they are known to have simpler waveforms compared to shallow earthquakes (e.g. Persh & Houston 2004), and the depth phases arrive sufficiently late so that they do not interfere with the duration measurements. In addition, Mw 8.3 Event 2 is of great interest, since it is the largest instrumentally-recorded deep earthquake to date (Ye et al.2013).

Figure 5.

Map view of the event locations and the rupture directions. The epicentral locations of Events 1 and 2, and the two major aftershocks of Event 2 shown with yellow, red and green GCMT mechanisms, respectively (Dziewoński et al.1981; Ekström et al.2012). The two aftershocks are both magnitude 6.7, and the southern event occurred on the same day (2013 May 24) as the main shock while the other happened on 2013 September 30. The yellow and red arrows show the horizontal projections of the inferred rupture extent and direction of Events 1 and 2, respectively, with three times amplification in length. The black contour lines with numbers are the slab depth contours of Slab 1.0 model (Hayes et al.2012). The background colour shows bathymetry and topography from ETOPO2 (distributed by National Centers for Environmental Information), and the inset map at the top left shows northeastern Pacific region with the plate boundaries (solid line; Bird 2003) and the dashed box corresponding to the area shown in the colour part of the figure.

Figure 5.

Map view of the event locations and the rupture directions. The epicentral locations of Events 1 and 2, and the two major aftershocks of Event 2 shown with yellow, red and green GCMT mechanisms, respectively (Dziewoński et al.1981; Ekström et al.2012). The two aftershocks are both magnitude 6.7, and the southern event occurred on the same day (2013 May 24) as the main shock while the other happened on 2013 September 30. The yellow and red arrows show the horizontal projections of the inferred rupture extent and direction of Events 1 and 2, respectively, with three times amplification in length. The black contour lines with numbers are the slab depth contours of Slab 1.0 model (Hayes et al.2012). The background colour shows bathymetry and topography from ETOPO2 (distributed by National Centers for Environmental Information), and the inset map at the top left shows northeastern Pacific region with the plate boundaries (solid line; Bird 2003) and the dashed box corresponding to the area shown in the colour part of the figure.

The measurements of duration are made using vertical broadband recordings available at the Data Management Center (DMC) of Incorporated Research Institution for Seismology (IRIS), and Full Range Seismograph Network (F-net) of Japan (Okada et al.2004). If a station has multiple vertical broadband sensors, one instrument is chosen so that one station location corresponds to one data point. Stations within 96° distance are selected to focus on the first-arriving compressional wave, that is, p or P phase, which results in 1187 stations for Event 1 and 1144 for Event 2. To avoid multiple arrivals within the first-arrival window, stations within the triplication distance range, that is, between 10$$\deg$$ and 20° based upon the IASP91 model (Kennett & Engdahl 1991), are removed (Fig. 6). Finally, stations with problematic instrument response, data segmentation, or low signal-to-noise ratio are discarded. After these considerations, the data set consists of 1102 and 1037 waveforms for Events 1 and 2, respectively (Fig. 6a). The upper-hemispheric coverage is poorer for Event 2 compared to Event 1, since most of the nearby stations in Japan lie within the triplication distance range. Before duration measurements are made, instrument responses are removed, and the seismograms are converted into displacements and bandpass filtered between 0.01 and 5 Hz using a four-pole two-pass Butterworth filter (Butterworth 1930).

Figure 6.

Station distribution and the triplication distance range. (a) Distribution of the stations used for the duration measurements of both Events 1 and 2 (black triangles), Event 1 only (yellow triangles) and Event 2 only (red triangles). The epicentres of the two events are shown by yellow (Event 1) and red (Event 2) stars. Yellow and red annuli represent the triplication distance ranges for Event 1 and Event 2, respectively. The map is plotted with azimuthal equidistant projection. (b) Dependence of ray parameter on distance for an event at 583-km depth calculated using a 1-D model IASP91 (Kennett & Engdahl 1991). The distance range of triplication is indicated by the shaded area. Solid and dashed curves correspond to P and p phases, respectively.

Figure 6.

Station distribution and the triplication distance range. (a) Distribution of the stations used for the duration measurements of both Events 1 and 2 (black triangles), Event 1 only (yellow triangles) and Event 2 only (red triangles). The epicentres of the two events are shown by yellow (Event 1) and red (Event 2) stars. Yellow and red annuli represent the triplication distance ranges for Event 1 and Event 2, respectively. The map is plotted with azimuthal equidistant projection. (b) Dependence of ray parameter on distance for an event at 583-km depth calculated using a 1-D model IASP91 (Kennett & Engdahl 1991). The distance range of triplication is indicated by the shaded area. Solid and dashed curves correspond to P and p phases, respectively.

The measured duration at a given station, τobs, is obtained using three time points, t1, t2 and t3 (Fig. 7a). The waveform onset time, t1, is relatively clear, but the end time can be difficult to choose. In order to incorporate this ambiguity, two end times, t2 and t3, are selected (Fig. 7a). Using these three picks, the duration is determined as

(21)
$$\tau ^{\text{obs}} = \frac{t_2+t_3}{2}-t_1,$$
and σ = (t3t2)/2 is assigned as its uncertainty. If the end time is clear so that t2 and t3 are identical, half of the shortest period for the filtered data, that is 0.1 s, is used as the uncertainty.

Figure 7.

A record section showing the directivity effect with respect to distance for Event 1. (a) Example displacement record from station KGM in southern Japan. The onset of the compressional wave is denoted by t1 and two values (t2 and t3) are selected as the end times. (b) Record section using stations with similar azimuthal angles (200°–250°, i.e. northern Japan and southeast Asia) aligned using hand-picked t1 times (grey line). The selected end times t2 and t3 are indicated by the grey swath.

Figure 7.

A record section showing the directivity effect with respect to distance for Event 1. (a) Example displacement record from station KGM in southern Japan. The onset of the compressional wave is denoted by t1 and two values (t2 and t3) are selected as the end times. (b) Record section using stations with similar azimuthal angles (200°–250°, i.e. northern Japan and southeast Asia) aligned using hand-picked t1 times (grey line). The selected end times t2 and t3 are indicated by the grey swath.

In addition to the measurement uncertainty, uneven data coverage needs to be considered. For both events, about 70 per cent of the stations are located within the United States, and this non-uniform distribution can lead to biased estimate of the rupture properties. Each datum is, therefore, weighted as

(22)
$$w = \frac{1}{N_{\rm d}\sqrt{\sigma } } ,$$
where σ is the uncertainty in the duration measurement and Nd is the number of stations within a given distance from the station. We choose $$\sqrt{\sigma }$$ instead of a more conventional σ2, because the inversion becomes heavily weighted toward only a few data with small σ if σ2 is used. The distances used to calculate Nd are determined on the focal sphere rather than the physical distances between stations, since coverage of the focal sphere is the one that matters (Fig. 8). In this study, distances of 3° and 5° are considered, and they yield similar results within the uncertainty bounds. Unless otherwise noted, the results described in the following sections are based upon the 3° bin size.

Figure 8.

Distribution of measured duration and uncertainty for Events 1 and 2. (a) Distribution of measured durations (coloured circles) plotted on the upper (top left) and lower (bottom left) hemispheres for Event 1. The uncertainties (right column) are shown in the same manner. Nodal planes (solid lines) are from the GCMT solution (Dziewoński et al.1981; Ekström et al.2012). The colour of each circle represents the measured duration or uncertainty. The arrow head and tail notation is used to show the inferred rupture direction. (b) Same as (a) except for Event 2. Note that the colour scale for the observed duration is different from (a) while the scheme for the uncertainty is the same. The grey shaded area indicates the range where take-off dip angle is 50° ± 15° which is used in Fig. 12(a).

Figure 8.

Distribution of measured duration and uncertainty for Events 1 and 2. (a) Distribution of measured durations (coloured circles) plotted on the upper (top left) and lower (bottom left) hemispheres for Event 1. The uncertainties (right column) are shown in the same manner. Nodal planes (solid lines) are from the GCMT solution (Dziewoński et al.1981; Ekström et al.2012). The colour of each circle represents the measured duration or uncertainty. The arrow head and tail notation is used to show the inferred rupture direction. (b) Same as (a) except for Event 2. Note that the colour scale for the observed duration is different from (a) while the scheme for the uncertainty is the same. The grey shaded area indicates the range where take-off dip angle is 50° ± 15° which is used in Fig. 12(a).

The measured durations projected onto the focal sphere show longest duration around the southwest part of the upper hemisphere and shortest duration around its antipode for Event 1. The gradual change in duration suggests that the unilateral rupture assumption is valid for this event. Event 2, on the other hand, shows more scatter than Event 1. There are regions, such as the southwestern quadrant, where long and short durations are both observed. This may partially be caused by the complex waveforms of the larger event and the difficulties in identifying the end time for Event 2 (Fig. 8). However, the longer durations for Event 2 are generally observed toward the northeast direction, whereas most of the short durations are seen on the southwest side of the lower hemisphere. Note that the duration distributions for Events 1 and 2 have reversed patterns (Fig. 8), suggesting that the ruptures propagated in opposite directions. Another feature that becomes obvious is the dip dependence. It can also be illustrated by a record section with respect to distance for a constant take-off azimuth (Fig. 7b), which supports the necessity of the 3-D analysis including the dip angle.

## RESULTS

### 2012 August 14, Mw 7.7 earthquake (Event 1)

Event 1 has observed durations ranging between 15 and 35 s with the average duration of 23 s from the 1102 stations. The misfit function obtained through the initial grid search is smoothly varying, and has a well-defined global minimum in dip and azimuth space with the preferred solutions of k0, γ0 and ϕ0 of 0.26, 30° and 40°, respectively (Fig. 9a). Combined with the 23-s average duration, this model gives a variance reduction of 99.3 per cent and χ2 of 7.6. The variance reduction is nearly 100 per cent since the average duration explains the main feature (finite duration) of the data. To check the fit to the fluctuation rather than the average value of the data, the variance reduction for the deviation of duration, that is, τi − τ0, is calculated. It gives 73.1 per cent for the starting model.

Figure 9.

Grid search and inversion results for Event 1. (a) Result of the initial grid search to determine the starting model. The grey shaded area is where k is negative and the white is where it is positive. The two minima in misfit are shown with crosses. (b) Iterative inversion result for the source duration τ (black) and the parameter k (red). Solutions of each parameter at each iteration stage are shown, where black circle and red square correspond to the parameters on the left and right axes, respectively. (c) Same as (b) except for the rupture dip γ and the rupture azimuth ϕ. (d) The misfit is calculated using eq. (11).

Figure 9.

Grid search and inversion results for Event 1. (a) Result of the initial grid search to determine the starting model. The grey shaded area is where k is negative and the white is where it is positive. The two minima in misfit are shown with crosses. (b) Iterative inversion result for the source duration τ (black) and the parameter k (red). Solutions of each parameter at each iteration stage are shown, where black circle and red square correspond to the parameters on the left and right axes, respectively. (c) Same as (b) except for the rupture dip γ and the rupture azimuth ϕ. (d) The misfit is calculated using eq. (11).

The next procedure, the iterative inversion, converges rapidly within three iterations (Figs 9b–d). In general, the misfit decreases most significantly at the first iteration, and becomes almost constant beyond the third iteration. The source duration reaches the final value with a single iteration, whereas the other parameters require three iterations to converge. The final solutions for the source duration, the parameter k, the dip and azimuth of rupture are 26 s, 0.27, 48° and 42°, respectively. The variance reduction and χ2 for this model are 99.5 per cent and 5.3, respectively, improving fits by 0.2 per cent and 2.3 compared to the starting model. The variance reduction for the duration fluctuation becomes 81.0 per cent, indicating better fit to the observed pattern in the data. The same set of final solutions are obtained even for extreme starting models, for example, where all four parameters are set to unity (i.e. source duration of 1 second, k of 1.0, and dip and azimuth of 1°). These tests with different starting models show that the final model does not depend upon the choice of the starting model. In order to assess the uncertainties of the solutions, bootstrapping is performed, where the data are sampled randomly and inverted to obtain a range of parameters that is warranted by the available data (Supporting Information). The bootstrapping analysis (Fig. S1) provides the preferred rupture model consisting of source duration of 26 ± 1 s, the k parameter of 0.27 ± 0.02, dip of 48° ± 8° and azimuth of 42° ± 5° (Table 1). The dip of the rupture is the least constrained of the four variables, which is expected from the limited dip coverage using the p and P phases (Fig. 8a).

Table 1.

Rupture parameters for the two deep earthquakes. Summary of the four parameters, that is, source duration (τ), k, rupture dip (γ) and azimuth (ϕ), as obtained from the iterative inversion. The rupture speed (VR) and extent (L) are calculated based upon the P-wave speed from the IASP91 model (Kennett & Engdahl 1991) and zero rise time assumption. The rupture speed is also given in terms of percentage of the shear-wave speed in parentheses. Events 1 and 2 correspond to an Mw 7.7 earthquake from 2012 August 14 and an Mw 8.3 earthquake from 2013 May 24, respectively. The same set of parameters are given for each episode of the two-episode analysis for Event 2, except that the timing (t) is shown instead of source duration. The episodes associated with clusters A and B of Event 2 are labelled as ‘2A’ and ‘2B’, respectively.

Event τ, t (s) k γ (deg) ϕ (deg) VR (km s−1L (km)
26 ± 1 0.27 ± 0.02 48 ± 8 42 ± 5 2.7 ± 0.2 (49 ± 4 per cent) 69 ± 6
37 ± 1 0.14 ± 0.02 − 19 ± 8 243 ± 11 1.4 ± 0.2 (26 ± 3 per cent) 52 ± 7
2A 33 ± 2 0.33 ± 0.10 − 22 ± 8 253 ± 9 3.3 ± 0.9 (61 ± 17 per cent) 111 ± 25
2B 33 ± 3 0.14 ± 0.10 − 15 ± 49 107 ± 58 1.4 ± 0.9 (26 ± 16 per cent) 47 ± 27
Event τ, t (s) k γ (deg) ϕ (deg) VR (km s−1L (km)
26 ± 1 0.27 ± 0.02 48 ± 8 42 ± 5 2.7 ± 0.2 (49 ± 4 per cent) 69 ± 6
37 ± 1 0.14 ± 0.02 − 19 ± 8 243 ± 11 1.4 ± 0.2 (26 ± 3 per cent) 52 ± 7
2A 33 ± 2 0.33 ± 0.10 − 22 ± 8 253 ± 9 3.3 ± 0.9 (61 ± 17 per cent) 111 ± 25
2B 33 ± 3 0.14 ± 0.10 − 15 ± 49 107 ± 58 1.4 ± 0.9 (26 ± 16 per cent) 47 ± 27

The constraints on the dip and azimuth indicate that the rupture propagates downdip towards the northeast. This direction is nearly parallel to the Kuril trench where the Pacific plate subducts underneath the North American plate (Fig. 5). Using the k parameter as the ratio of rupture to compressional wave speed and the 9.9 km s−1 compressional wave speed from the IASP91 model at the hypocentral depth of 583 km (Kennett & Engdahl 1991), the rupture speed is inferred to be 2.7 ± 0.2 km s−1 (Table 1). This propagation speed corresponds to 49 ± 4 per cent of the shear wave speed, using 5.4 km s−1 for shear wave speed at the hypocentral depth (Table 1). Note that these speed estimates are obtained under the assumption of zero rise time, and should be considered as the lower bound. Combining the source duration of 25 ± 1 s and the rupture speed estimate, the extent of this event is 69 ± 6 km (Table 1). These results do not change significantly if the 5° distance weighting scheme is used, and the preferred values from the larger bin fall within the uncertainties of the 3°-bin results.

To complement the directivity analysis, we perform back-projection modelling of the North American data obtained from IRIS DMC with 12-station supplementation from the Superior Rifting EarthScope Experiment (Stein et al.2011; Supporting Information). Due to limited spatial and temporal resolution of the back-projection technique, tracking rupture for spatially compact events such as the one studied here becomes difficult (Kiser & Ishii 2012); hence we focus on identifying subevents using the coherency function analysis (Ishii 2011; Supporting Information). The back-projection analysis identifies three subevents (Fig. 10; Table S1). The first subevent is imposed to be consistent with the hypocentral location while the second subevent is located to the east of the first event, followed by a third subevent which is to the north of the second subevent (northeast of the epicentre; Fig. 10a). The main high-frequency energy release occurs during the second and third subevents, and the total duration for the earthquake is about 20 s (Fig. 10b). If we assume that the rupture propagated smoothly from one subevent to another, the rupture speed between subevents 1 and 2 is about 2.7 km s−1 and between subevents 2 and 3 is about 3.2 km s−1. Taking the location and timing of the first and last subevents, the ‘average’ rupture speed is about 2.1 km s−1. Note that these estimates assume propagation on a horizontal plane, and as such, represent lower limits.

Figure 10.

Back-projection results of the Mw 7.7 2012 August 14 earthquake. (a) Summary of the spatial variation. The background colour shows the time-integrated coherency function over the duration of the event where warm colours indicate high coherency and cold colours indicate low coherency. The black contours are plotted at 10 per cent interval starting at 50 per cent of the highest coherency value. The large white star shows the epicentral location as given by the NEIC catalogue. There are three coloured stars with black outline corresponding to the three subevents identified (Supporting Information Table S1). The red star which is hidden behind the white epicentral star is at the beginning of the event, the yellow star is the middle subevent and the blue star indicates the last subevent location. The white arrow shows the average direction towards the stations used in the analysis. (b) Summary of the temporal variation with the black curve showing the maximum coherency value as a function of time. The coloured vertical bars correspond to the times at which the three coloured subevents are shown in (a). The grey curve shows the normalized linear-stacking maximum values as a function of time.

Figure 10.

Back-projection results of the Mw 7.7 2012 August 14 earthquake. (a) Summary of the spatial variation. The background colour shows the time-integrated coherency function over the duration of the event where warm colours indicate high coherency and cold colours indicate low coherency. The black contours are plotted at 10 per cent interval starting at 50 per cent of the highest coherency value. The large white star shows the epicentral location as given by the NEIC catalogue. There are three coloured stars with black outline corresponding to the three subevents identified (Supporting Information Table S1). The red star which is hidden behind the white epicentral star is at the beginning of the event, the yellow star is the middle subevent and the blue star indicates the last subevent location. The white arrow shows the average direction towards the stations used in the analysis. (b) Summary of the temporal variation with the black curve showing the maximum coherency value as a function of time. The coloured vertical bars correspond to the times at which the three coloured subevents are shown in (a). The grey curve shows the normalized linear-stacking maximum values as a function of time.

Given that the inversion results from the directivity analysis are in 3-D while the back-projection results are in 2-D, the parameters that can be directly compared without horizontal projection are source duration and azimuth. The source duration from back-projection is about 6 s shorter than the preferred value from the directivity analysis. This may be caused by the difference between 2-D and 3-D analyses, or by the difference in the frequency of the data used for the two analyses—for example, back-projection may not be capturing the low-frequency energy near the end of the rupture. On the other hand, the azimuthal direction obtained from the inversion is consistent with the back-projected location of the last subevent, that is, subevent 3, with respect to the epicentre. The azimuth of subevent 3 is 44° (northeast) from the epicentre, and this angle is within the uncertainty of the preferred azimuth from the directivity analysis (42° ± 5°).

When the rupture extent and speed from the inversion are projected onto a horizontal plane using the estimated dip, they can be compared with the back-projection results. During this process, the large uncertainty in the dip angle propagates into the projected estimates, and as a result, the rupture speed and extent acquire large uncertainties, 1.8 ± 0.3 km s−1 and 46 ± 8 km, respectively. The horizontal propagation speed is consistent with the ‘average’ rupture speed from back-projection of 2.1 km s−1. The horizontal rupture extent also compares well to the distance between the epicentre and the last subevent from the back-projection analysis, although the former is about 13 km longer. This difference arises from longer source duration obtained from the directivity analysis, since the rupture extent is based upon the source duration.

Considering that Event 1 shows a significant rupture dip, the event can be used to verify the importance of the dip in the directivity analysis that is emphasized in Section 2.2. If the inversion is performed with the dip fixed to zero, it yields 22 s, 0.26 and 32° for the source duration, parameter k, and rupture azimuth, respectively. Although the value of k remains similar to the original value of 0.27, the rupture azimuth is about 10° smaller compared to the preferred solution. The source duration is also underestimated, and becomes closer to the back-projected value. The amount of underestimation of 3 s is consistent with our prediction for the downdip rupture case (Fig. 3b) and the bias calculated using eq. (9).

### 2013 May 24, Mw 8.3 earthquake (Event 2)

The magnitude 7.7 2012 August 14 (Event 1) is followed 9 months later by Event 2, which is 0.6 magnitude units larger. As the largest deep event ever recorded, Event 2 has longer observed durations compared to Event 1, ranging between 27 and 55 s with large uncertainties (Fig. 8). The measurements also contain significant scatter, suggesting that the unilateral rupture assumption may not be valid. Nonetheless, we first perform the inversion to examine how well the data can be modelled with the unilateral assumption. We then apply the complex rupture model to investigate whether multiple subevents can be identified.

Noting that the observed durations are not smoothly varying, have large uncertainties, and the upper-hemispheric coverage is poor, the inversion scheme is modified; instead of using the average duration as the initial source duration, it is included as an additional parameter in the grid search. The source duration is varied between the minimum and maximum observed values with 5-s increment. The solutions obtained through this modified grid search are 35 s for the source duration (in contrast to the average duration of 38 s), 0.18 for the k parameter, −40° dip and 240° azimuth. This starting model gives variance reduction of 99.2 per cent and χ2 of 17.2. The variance reduction that measures the fit to the data fluctuation, as discussed for Event 1, is 30.0 per cent.

The subsequent iterative inversion using this starting model converges after three iterations, similar to the behaviour observed for Event 1. The largest perturbation occurs in the first iteration, and the model parameters become nearly constant beyond the third iteration. The model values after three iterations are 37 s, 0.14, −19° and 243° for the source duration, k, dip, and azimuth, respectively. However, the fit to data does not change substantially with the iterative inversion. The total change is a 0.03 per cent increase in the variance reduction and 0.6 decrease in the χ2. The variance reduction for the duration fluctuation becomes 32 per cent, a 2 per cent improvement compared to the starting model.

The results of the iterative inversion for this event do not depend strongly on the initial model. Starting with different models produces identical results, but convergence is slower when the initial model is far from the final solution. If all the parameters in the starting model are set to unity, it takes more than ten iterations to converge. The slow convergence is mainly caused by difficulties in reaching a reasonable value for the source duration using the perturbation approach. If the iterative inversion starts with a source duration that is relatively close to the final solution, the inversion converges after about five iterations regardless of the initial values of k, dip, and azimuth. These tests indicate that Event 2 is more sensitive than Event 1 to the starting model parameters, especially to the initial source duration, although the final solutions are nonetheless robust.

The estimates of uncertainty based upon the bootstrapping procedure are also influenced by the scatter in the data and poor coverage (Supporting Information; Fig. S2). Recognizing that the two measurements in the upper hemisphere are an important part of the data set (Supporting Information), the bootstrapping procedure is carried out with the two upper-hemispheric data included in every set. This results in nearly Gaussian distributions for all four parameters, and the preferred solutions are 37 ± 1 s, 0.14 ± 0.02, −19° ± 8° and 243° ± 11° for the source duration, k parameter, the dip and the azimuth, respectively (Table 1). The rupture speed of this event, calculated with zero rise time assumption and the compressional wave speed given in IASP91 at 602-km depth (Kennett & Engdahl 1991), is 1.4 ± 0.2 km s−1, corresponding to about 26 ± 3 per cent of shear wave speed. Combining this estimate with the source duration gives the rupture extent of 52 ± 7 km (Table 1). These results are weakly dependent upon the size of the distance bin used to define the weighting values. The dip parameter changes to −31° ± 5° if the 5° weighting scheme is used, which is a consequence of the poor upper-hemispheric coverage. However, the 5°-bin results for the other parameters fall within the uncertainty of the 3°-bin model.

For this earthquake, the back-projection analysis shows a complex rupture pattern with six subevents (Fig. 11; Table S2). The initial location, subevent 1, is imposed to be consistent with the hypocentral location. Subevents 2 and 3 occur close in time, separated by about 2 s, but at significantly different locations, northeast and south of the hypocentre, respectively. Their locations, combined with snapshot images, suggest that the initial rupture from the epicentre is bilateral. Assuming that the rupture is continuous between subevents 1 and 2, and between subevents 1 and 3, the corresponding rupture speeds are 3.6 km s−1 and 2.3 km s−1, respectively. Subsequent subevents 4, 5 and 6 appear isolated, and rupture propagation between subevents (e.g. rupture from subevent 3 to 4) is difficult to resolve. Subevents 4 and 5 are located to the southeast of the epicentre and are associated with the largest high-frequency energy release (Fig. 11). The last subevent, subevent 6, is to the southwest of the epicentre. The duration of the entire event is estimated to be about 38 s (Supporting Information).

Figure 11.

Back-projection results of the Mw 8.3 2013 May 24 earthquake. (a) Same as in Fig. 10(a) except that there are six subevent stars (Supporting Information Table S2). The order of the subevent occurrence is red, orange (northeast of the epicentre), yellow (south of the epicentre), green, dark green and blue (southwest of the epicentre). The grey circles show the aftershocks available in the NEIC catalogue between 2013 May 24 and 2014 June 1. Lighter grey indicates shallower depth and darker grey indicates deeper depth with minimum and maximum depths of 509 and 642 km, respectively. (b) Same as in Fig. 10(b). For resolvability of subevent 3 (yellow vertical bar), see Supporting Information.

Figure 11.

Back-projection results of the Mw 8.3 2013 May 24 earthquake. (a) Same as in Fig. 10(a) except that there are six subevent stars (Supporting Information Table S2). The order of the subevent occurrence is red, orange (northeast of the epicentre), yellow (south of the epicentre), green, dark green and blue (southwest of the epicentre). The grey circles show the aftershocks available in the NEIC catalogue between 2013 May 24 and 2014 June 1. Lighter grey indicates shallower depth and darker grey indicates deeper depth with minimum and maximum depths of 509 and 642 km, respectively. (b) Same as in Fig. 10(b). For resolvability of subevent 3 (yellow vertical bar), see Supporting Information.

The source durations obtained from both directivity and back-projection analyses are in good agreement. The rupture extent and speed from the directivity analysis, projected onto a horizontal plane using the dip estimate, become 49 ± 7 km and 1.3 ± 0.2 km s−1, respectively. These are compatible with 50-km distance between the epicentre and subevent 6, and the ‘average’ rupture speed of 1.3 km s−1 calculated using the distance and duration. As in the case of Event 1, our rupture azimuth for Event 2 (243° ± 11°) is consistent with the location of the last subevent (234°) from the epicentre.

Noting that the back-projection analysis indicates complex rupture pattern, we apply multi-episode directivity analysis to Event 2. To demonstrate that this approach is warranted, the duration measurements with take-off angles of about 50° (Fig. 8b) are binned and plotted as a function of azimuth (Fig. 12a). While the values are mostly consistent with the observed timing from our inversion results, the dependence on azimuth is not a smoothly-varying sinusoid but has a cusp around 180°, suggesting that at least two episodes are involved. Assuming that there are two episodes contributing to the observed end timing of the P waveform, the simulated annealing algorithm (e.g. Metropolis et al.1953; Kirkpatrick et al.1983; Ingber 1995) is applied. For each episode, there are four model parameters (episode timing, k, dip and azimuth) as in the unilateral directivity analysis, resulting in eight model parameters being sought. The parameter k for the nth episode is expressed as Ln/(tnV), where Ln is the distance between the episode and the hypocentre, tn is the timing with respect to the hypocentral time and V is the seismic wave speed at the source (eq. 19). The simulated annealing algorithm requires initial input for a variable, temperature, associated with the algorithm's search mechanism, and it is set to 500. This is sufficient to sample the model space to identify the globally optimal solutions, but other values between 100 and 1000 produce similar outcomes. Simulated annealing is performed one thousand times starting from models that are randomly selected from the possible range of values: between 27 and 55 s for episode timings, 0 and 1 for k, −90° and 90° for dip angles, and 0° and 360° for azimuths.

Figure 12.

Results of the two-episode analysis of the Mw 8.3 2013 May 24 earthquake. (a) Measured duration (grey dots) as a function of azimuth for those that have take-off dip angle of 50° ± 15°, that is, data points within the grey annulus shown in Fig. 8(b). White circles with black outline indicate data that have been binned for every 30° in azimuth. These data are compared to predicted durations from the 284 solutions (pink curves) and the best solution (purple curve). The prediction from the unilateral analysis is shown by the blue curve. (b) Epicentral locations of the 284 solutions from the simulated annealing. Clusters A (circles) and B (crosses) are coloured by the relative timing (tAtB) of the pair of episodes, and the sizes of the symbols are proportional to the fit to the data with larger symbols indicating better fit. Purple circle and cross show the position of the best solution. The NEIC hypocentral location for the event is shown by the star. (c) Same as in (b) except that the solutions are shown in the depth–longitude space. (d) Two-episode timing (top), distance (second from top), dip (second from bottom) and azimuth (bottom) values obtained from the simulated annealing as functions of the timing between the two episodes. The symbols are the same as in (b) except that cluster A is further emphasized by black colour and cluster B by grey colour.

Figure 12.

Results of the two-episode analysis of the Mw 8.3 2013 May 24 earthquake. (a) Measured duration (grey dots) as a function of azimuth for those that have take-off dip angle of 50° ± 15°, that is, data points within the grey annulus shown in Fig. 8(b). White circles with black outline indicate data that have been binned for every 30° in azimuth. These data are compared to predicted durations from the 284 solutions (pink curves) and the best solution (purple curve). The prediction from the unilateral analysis is shown by the blue curve. (b) Epicentral locations of the 284 solutions from the simulated annealing. Clusters A (circles) and B (crosses) are coloured by the relative timing (tAtB) of the pair of episodes, and the sizes of the symbols are proportional to the fit to the data with larger symbols indicating better fit. Purple circle and cross show the position of the best solution. The NEIC hypocentral location for the event is shown by the star. (c) Same as in (b) except that the solutions are shown in the depth–longitude space. (d) Two-episode timing (top), distance (second from top), dip (second from bottom) and azimuth (bottom) values obtained from the simulated annealing as functions of the timing between the two episodes. The symbols are the same as in (b) except that cluster A is further emphasized by black colour and cluster B by grey colour.

The best solution from the two-episode analysis consists of episodes that occur almost at the same time, 33 and 34 s from the hypocentral time. The parameter k, distance, dip, and azimuth obtained are 0.33, 112 km, −21° and 257°, respectively, for one of the episodes, which we refer to as ‘episode A’, and 0.10, 33 km, −20° and 114°, respectively, for the other, ‘episode B’ (Figs 12b and c). This solution gives χ2 of 13.8 and the variance reduction for the duration fluctuation of 45 per cent. The χ2 and the variance reduction are 16.6 and 32 per cent, respectively, for the model based upon unilateral assumption; hence the improvement in fit using the two-episode model is statistically significant (Fig. 12a). Among the thousand simulated annealing solutions, 917 fit the data statistically better than the unilateral model, that is, the χ2 value that incorporates the increase in the number of model parameters improves (Fig. 12a). Of the 917 solutions, 633 have variance reduction for the duration fluctuation that is above 41 per cent, that is, within 90 per cent of the variance reduction obtained for the best solution of which 284 are within 95 per cent. Not surprisingly, these solutions give similar values for model parameters, and the 3-D locations form two distinct clusters (Figs 12b and c). Here, we focus on the 284 solutions to examine the distribution and range of model parameters.

One cluster that includes episode A is located to the southwest of the hypocentre, and this set of solutions is hereafter referred to as ‘cluster A’. It has dip and azimuth angles ranging from −37° to −12° and from 243° to 265°, respectively (Figs 12b–d). They are consistent with the solution from the unilateral analysis, as well as the azimuth of the final subevent from the back-projection analysis (Fig. 11a). The timing and distance of solutions in this cluster range between 30 and 35 s and between 83 and 150 km, respectively. Compared to the results from the unilateral analysis, the timing is 2–6 s earlier and the distance is about 30–100 km longer. This cluster of solutions reproduce durations observed at stations located at about 0°–180° azimuth (Fig. 12a). The other cluster around episode B of the best solution, referred to as ‘cluster B’, is located east of the epicentre. It is a new feature not included in the previous unilateral rupture analysis. Its timing varies between 30 and 37 s, and its location distribution spreads over a larger volume than that of cluster A (Figs 12b and c). The distance from the hypocentre, dip, and azimuth range from 16 to 85 km, from −66° to 48° and from 51° to 194°, respectively (Fig. 12d). By introducing cluster B, the observed data pattern between about 180° and 330° in station azimuth can be explained, and the combination of clusters A and B generates the cusp around 180° (Fig. 12a).

There is a trade-off between timing of episodes A and B (Fig. 12d). A solution in cluster A that occurs early has its counterpart in cluster B that occurs late, and vice versa. The relative timing (tAtB) varies from −6 to 4 s, implying that the precedence of one episode over another is not constrained with currently available data. The relative timing also influences the 3-D locations of the two episodes. As the relative timing increases, the solutions in cluster A move toward southeast and downdip, that is, they move closer to the hypocentre. In contrast, episode locations of cluster B change toward northeast and updip, that is, they move farther away from the hypocentre (Figs 12b and c). These trade-offs complicate comparison with back-projection results. The cluster A is located in the direction of the last subevent identified through back-projection analysis. The cluster B, however, is not compatible with the penultimate subevent, subevent 5 (Fig. 11). The episodes within cluster B are more consistent with the location of subevent 4 that occurs 15 s before the last subevent. This discrepancy can be reconciled if there is significant low-frequency energy release near the subevent 4 location that lasts for some time (Xu et al.2009; Kiser & Ishii 2011).

As in the unilateral directivity analysis, our preferred two-episode model for Event 2 is obtained using the distributions of parameters from the best 284 solutions, which are well-approximated as Gaussian. The timing, k, distance, dip, and azimuth are 33 ± 2 s, 0.33 ± 0.10, 111 ± 25 km, −22° ± 8°, 253° ± 9°, respectively, for the episode associated with cluster A, and 33 ± 3 s, 0.14 ± 0.10, 47 ± 27 km, −15° ± 49°, 107° ± 58°, respectively, for the other episode corresponding to cluster B (Table 1). In general, uncertainties for cluster B are larger than those of cluster A. This is caused by the fact that the amount of data used to constrain cluster B is less than that for cluster A—that is, a large fraction of the data is explained by cluster A (which also explains why the unilateral model is compatible with cluster A). The total source duration of Event 2 is calculated for each of the 284 solutions as the maximum of the times of the two episodes, and the distribution of the values results in the estimate of 34 ± 2 s (Table 1). It is similar but shorter than 37 ± 1 s from the unilateral model, consistent with the fact that unilateral assumption for a complex rupture overestimates the source duration. The episode timings and locations can be used to calculate corresponding rupture speeds, and they are 3.3 ± 0.9 km s−1 for cluster A and 1.4 ± 0.9 km s−1 for cluster B (Table 1). Note that these estimates assume zero rise time, and are associated with distances from the hypocentre to each episode, and thus, are ‘averaged’ speeds, and may not reflect the true speeds.

To test if additional episodes are required to improve fit to the data, we apply three-episode analysis. The results show no significant improvement in the misfit, and the locations of the third episode obtained from different simulated annealing solutions do not form a cluster. Hence, including three episodes in the analysis does not help to explain the data more than the two-episode analysis does.

## DISCUSSION

Deep earthquakes are known for their lack of aftershocks (e.g. Frohlich 1989), but the Mw 8.3 earthquake from 2013 May 24 (Event 2) is followed by a sequence of aftershocks including two at magnitude 6.7 (Figs 5 and 11a). Although previous studies have suggested that aftershock locations of deep main shocks encompass a diffuse volume rather than a single plane (Tibi et al.2003; Chen et al.2014), the aftershock distribution of Event 2 is nearly planar (Fig. 11a). The azimuths of the aftershocks with respect to the mainshock hypocentre are compatible with the southern direction resolved as subevents 3, 4 and 5 in back-projection analysis, or as a subset of cluster B in the two-episode analysis. The depths of the aftershocks are nearly constant, indicating a plane that is subhorizontal rather than subvertical, consistent with the directivity solutions.

For both Events 1 and 2, our inverted source durations (26 ± 1 and 37 ± 1 s, respectively) are significantly shorter than those from the Global Centroid Moment Tensor (GCMT) solutions (35.6 and 71.4 s; Dziewoński et al.1981; Ekström et al.2012). The discrepancy can be attributed to the magnitude-duration relationship used to calculate the source durations for GCMT solutions (Ekström & Engdahl 1989; Ekström et al.2012), that is,

(23)
$$\tau _{{\rm h}} = 1.05 \times 10^{-8} M_0^{1/3},$$
where τh is the half duration in seconds and M0 is the scalar moment in dyne cm. This magnitude-duration relationship has been empirically determined using events at shallow depths (Ekström & Engdahl 1989), and there have been indications that deep-focus earthquakes deviate from this relationship and have shorter durations (e.g. Vidale & Houston 1993; Bos et al.1998; Persh & Houston 2004). The SCARDEC solution obtained using relatively high-frequency data (1–3 Hz; Vallée et al.2011) estimates the source time function of about 21 s for Event 1, and 38 s for Event 2. Event 2 has also been studied by other research groups, and durations between 30 and 40 s have been reported (Wei et al.2013; Ye et al.2013; Chen et al.2014; Kennett et al.2014; Zhan et al.2014). The durations obtained from the directivity analyses are more compatible with these estimates than the GCMT durations.

Another parameter with implications for the rupture mechanism is the propagation speed (e.g. Frohlich 1989). Our models based upon unilateral rupture assumption give rupture speeds of 2.7 ± 0.2 and 1.4 ± 0.2 km s−1, or 49 ± 4 and 26 ± 3 per cent of the shear wave speed, for Events 1 and 2, respectively. Deep-focus earthquakes are observed to have rupture propagation at 30 to 90 per cent of the shear wave speed (e.g. Frohlich 2006; Suzuki & Yagi 2011), and the estimate for Event 1 is well within this range, while the value for Event 2 is close to the lower bound. This does not, however, reflect peculiarities in the rupture process, but an inappropriate assumption in the modelling approach. Previous studies (Wei et al.2013; Ye et al.2013; Chen et al.2014; Meng et al.2014; Zhan et al.2014) and the back-projection analysis presented in this paper show that the Mw 8.3 earthquake ruptured in a complex pattern, which implies that the duration is longer than it would have taken from the hypocentre to the last part of the energy release. The directivity analysis only provides the total duration and last episode location with respect to the hypocentre; hence it will inevitably result in underestimated rupture speed. The back-projection analysis for this event shows rupture speeds of 2.3 and 3.6 km s−1 (on the horizontal plane) for initial parts of the rupture process, in good agreement with values suggested for deep earthquakes (e.g. Frohlich 2006; Suzuki & Yagi 2011).

Other parameters obtained in this study are associated with the geometry of the rupture plane, that is, the dip and azimuth. The rupture direction of Event 1 is inferred to be dipping by 48° ± 8° with azimuth of 42° ± 5° while that of Event 2 is nearly horizontal with almost opposite azimuth. When compared to the nodal planes of the GCMT solution (Dziewoński et al.1981; Ekström et al.2012), the rupture direction obtained for Event 1 agrees with the nodal plane striking at 243° and dipping at 63°. For Event 2, GCMT has two nodal planes with 189° strike and 11° dip, and 12° strike and 79° dip. Dip angle estimates from the unilateral and two-episode analyses of this study favour the first nodal plane. The preference for a nearly horizontal plane is consistent with the result of Wei et al. (2013) where the authors examine the fit to data using the two nodal planes of the NEIC W-phase focal mechanism. Meng et al. (2014) and Zhan et al. (2014) also argue for the horizontal plane based upon their back-projection analysis and subevent locations, respectively. Similarly, Ye et al. (2013) prefers the shallowly dipping plane although their data do not preclude a nearly vertical plane. On the other hand, Chen et al. (2014) suggests significant vertical extent using the difference in the durations of P and pP. However, as the authors point out, the onset and end times of pP phase are not as clear as those for the P phase, and the duration measurements, and thus the vertical rupture extent, may have large uncertainties. There is another study showing vertical trajectory of high-frequency energy propagating downward (Kennett et al.2014), and these authors propose that the discrepancy between their result and other studies may come from the difference in locations of the high- and low-frequency energy release.

Although the dip angle for the Mw 8.3 event is not as well-constrained as for the Mw 7.7 event, the difference is significant. Event 2 ruptures on a nearly horizontal plane (−19° ± 8°) while Event 1 ruptures with substantial dip (48° ± 8°). This change in geometry within the same subducting slab, separated only by 800 km, suggests the existence of a complex stress field. Seismic tomography studies of this subduction zone have shown complicated slab structure in this region. Some portions of the slab are observed to flatten around 650-km depth (e.g. Fukao et al.1992; van der Hilst et al.1993; Takenaka et al.1999) while other parts appear to penetrate into the lower mantle (e.g. Creager & Jordan 1984; Zhou et al.1990; Gaherty et al.1991; Schwartz et al.1991; Glennon & Chen 1993; Ding & Grand 1994; Pankow & Lay 1999). The differences in the dip angle between the two mainshocks considered in this study may reflect differences in the slab dynamics.

For the azimuthal constraint on Event 2, previous studies conclude that the dominant rupture direction is south or southeast, which is about 60° to 80° off from 243° ± 11° of the model presented in this manuscript. This discrepancy may be caused by the fact that the last subevent has small amplitude (Fig. 11; Supporting Information) that may have been obscured in other back-projection analyses. In our directivity analysis, the end time measurements emphasize the final slip location and the method is less prone to amplitude variations. Furthermore, south or southeast locations are included in the cluster of solutions for the second source in the two-episode analysis, compatible with models from other studies. The two-episode model also helps explain an observation that a line cannot fit the directivity pattern assuming a single rupture direction (e.g. Zhan et al.2014). In contrast, Chen et al. (2014) consider durations of P phase to obtain an azimuth of about 160° by finding the shortest duration. As described in Section 2.4, for a complicated rupture such as Event 2, the shortest observed duration is not the rupture propagation direction (Fig. 4). For these complex events, it is more effective to search for maxima in observed durations that correspond to the antipodal directions from different episodes.

## CONCLUSIONS

We have introduced and implemented a technique to obtain some earthquake source parameters based upon the 3-D directivity formulation, emphasizing the importance of dip angle in the analysis. For example, relying upon teleseismic data, that is, sparse sampling of the upper focal hemisphere, can result in a source duration that is about 20 per cent under or overestimated. Using the unilateral rupture model, our method requires no additional assumptions or a priori constraints such as fault plane orientation. The inversion is nonlinear, but a reasonable starting model can be obtained easily with a grid search, and the iterative part of inversion converges rapidly.

The inversion scheme seeks four rupture parameters—the source duration, parameter k which is related to rupture extent and speed, dip, and azimuth of propagation—and has been applied to two deep-focus events in the Sea of Okhotsk region, an Mw 7.7 event that occurred on 2012 August 14, and an Mw 8.3 event from 2013 May 24. The measured compressional-wave durations of the first event show smooth dependence on the take-off angle and azimuth, consistent with the expected pattern from a unilateral rupture propagation. The resulting rupture parameters are source duration of 26 ± 1 s, rupture speed of at least 2.7 ± 0.2 km s−1 (i.e. 49 ± 4 per cent of the shear wave speed), rupture extent of 69 ± 6 km, and propagation direction of 48° ± 8° dip and 42° ± 5° azimuth, nearly parallel to the Kuril trench. The second, larger earthquake, on the other hand, shows considerable complexity in the data, suggestive of a rupture that is not unilateral. The preferred solution has a source duration of 37 ± 1 s, rupture speed of at least 1.4 ± 0.2 km s−1 (i.e. 26 ± 3 per cent of the shear wave speed), and rupture direction of −19° ± 8° dip and 243° ± 11° azimuth. The rupture propagation is parallel to the Kuril trench, but is in the opposite direction to that of the Mw 7.7 event. We compare these results with back-projection analysis of P-wave data recorded in North America. The images show complex rupture behaviour for the second earthquake with at least three different directions of propagation, and that the directivity analysis with the unilateral rupture assumption gives values associated with the location of the last subevent, or the last part of the energy release.

In view of the complexities observed for the Mw 8.3 earthquake, we extend our approach beyond unilateral rupture by replacing the duration equation with a general expression for rupture with multiple episodes in arbitrary directions. This formulation provides insight into some challenges associated with directivity analysis. For example, the shortest observed duration is no longer the direction of rupture, and the assumption of unilateral rupture for a complex earthquake leads to overestimation of the source duration. Application of this method using the simulated annealing algorithm results in a two-episode solution that significantly improves the fit to the data. One cluster of locations is compatible with the location inferred from the unilateral modelling as well as the final subevent identified through back-projection analysis, and the other group of locations is close to some subevent locations from the back-projection analysis. The source duration estimate obtained from the multi-episode method is about 3 s shorter, consistent with the inherent overestimation of the source duration in the unilateral rupture model.

Improvements in the estimates can be obtained with stations that are more uniformly distributed in the take-off angle and azimuth. Achieving better coverage through installation of more seismic stations around the world is difficult, but the coverage, especially in the upper focal hemisphere, can be improved by incorporating other seismic phases such as the depth phase. For a complex earthquake involving multiple subevents, if one can identify end or peak-amplitude times of each subevent, they can be used to constrain subevent timing and location, as well as the overall rupture area. Some studies have followed this approach, and inferred the subevent locations (e.g. Bezada & Humphreys 2012; Wei et al.2013), but it is important to note that arrivals associated with two or more subevents do not appear as simple sinusoidal functions. Other efforts to automate directivity analysis without picking arrival or end times using the waveform stretching and matching approach (Warren & Silver 2006) can be combined with our inversion method. For complex non-unilateral ruptures, different amplitude and stretching associated with various subevents must and can be taken into account. Finally, the inversion scheme described in this paper can be applied to smaller events, possibly as small as Mw 5.7 (Warren & Shearer 2006), or shallower events as long as variations in body-wave durations are observed without interference from other phases such as pP.

The authors thank Yehuda Ben-Zion and two anonymous reviewers for comments and suggestions that helped improve the manuscript. We also thank Arizona Geological Survey, NOAA/NWS/West Coast-Alaska Tsunami Warning Center, Geoscience Australia, Alaska Volcano Observatory, IGPP (University of California, San Diego), Swiss Seismological Service, California Institute of Technology, Geological Survey of Canada, Geophysical Institute Academy of Science of the Czech Republic, Dublin Institute for Advanced Studies, Institut de Physique du Globe de Paris, GFZ Potsdam, Hong Kong Observatory, Central-Asian Institute of Applied Geosciences, Institute of Seismology of National Academy of Sciences of Kyrgyz Republic, Kazakhstan National Data Center, Lamont-Doherty Earth Observatory of Columbia University, Montana Bureau of Mines and Geology, Istituto Nazionale di Geofisica e Vulcanologia (Rome, Italy), Institut de Recherche pour le Developpement/Research Institute for Development (New Caledonia), Weston Observatory of Boston College, Saint Louis University Earthquake Center, NORSAR (Norway), Instituto Nicaragüense de Estudios Territoriales, Central Institute for Meteorology and Geodynamics (Austria), Plate Boundary Observatory (EarthScope), Institute of Geophysics (Polish Academy of Sciences), Regional Integrated Multi-Hazard Early Warning System for Africa and Asia, Romanian National Data Center, New Mexico Tech, Institute of Earth Sciences (Academia Sinica, Taiwan), University of Oregon, University of Utah, University of Washington, National Research Institute for Earth Science and Disaster Prevention in Japan, U.S. Geological Survey, National Science Foundation, the Data Management Center of Incorporated Research Institutions for Seismology, and S. van der Lee for providing the data. Some figures have been generated using the Generic Mapping Tools (Wessel & Smith 1991). SP is supported by the Samsung Scholarship and the Peirce fellowship.

## REFERENCES

Ammon
C.J.
et al
2005
Rupture process of the 2004 Sumatra-Andaman earthquake
Science

308
1133
1139
Beck
S.L.
Silver
P.
Wallace
T.C.
James
D.
1995
Directivity analysis of the deep Bolivian earthquake of June 9, 1994
Geophys. Res. Lett.

22
2257
2260
M.J.
Humphreys
E.D.
2012
Contrasting rupture processes during the April 11, 2010 deep-focus earthquake beneath Granada, Spain
Earth planet. Sci. Lett.

353
38
46
Bird
P.
2003
An updated digital model of plate boundaries
Geochem. Geophys. Geosyst.

4
1027
doi:10.1029/2001GC000252
Bos
A.G.
Nolet
G.
Rubin
A.
Houston
H.
Vidale
J.E.
1998
Duration of deep earthquakes determined by stacking of Global Seismograph Network seismograms
J. geophys. Res.

103
21 059
21 065
Butler
R.
Stewart
G. S.
Kanamori
H.
1979
The July 27, 1976 Tangshan, China earthquake–a complex sequence of intraplate events
Bull. seism. Soc. Am.

69
207
220
Butterworth
S.
1930
On the theory of filter amplifiers
Wireless Engineer

7
536
541
Caldeira
B.
Bezzeghoud
M.
Borges
J.F.
2010
DIRDOP: a directivity approach to determining the seismic rupture velocity vector
J. seismol

14
565
600
Chen
Y.
Wen
L.
Ji
C.
2014
A cascading failure during the 24 May 2013 great Okhotsk deep earthquake
J. geophys. Res.

119
3035
3049
Convertito
V.
Caccavale
M.
De Matteis
R.
Emolo
A.
Wald
D.
Zollo
A.
2012
Fault extent estimation for near-real-time ground-shaking map computation purposes
Bull. seism. Soc. Am.

102
661
679
Creager
K.C.
Jordan
T.H.
1984
Slab penetration into the lower mantle
J. geophys. Res.

89
3031
3049
Ding
X.-Y.
Grand
S.P.
1994
Seismic structure of the deep Kurile subduction zone
J. geophys. Res.

99
23 767
23 786
Dziewoński
A.M.
Chou
T.A.
Woodhouse
J.H.
1981
Determination of earthquake source parameters from waveform data for studies of global and regional seismicity
J. geophys. Res.

86
2825
2852
Ekström
G.
Engdahl
E.
1989
Earthquake source parameters and stress distribution in the Adak Island region of the central Aleutian Islands, Alaska
J. geophys. Res.

94
15 499
15 519
Ekström
G.
Nettles
M.
Dziewoński
A.M.
2012
The global CMT project 2004–2010: centroid-moment tensors for 13,017 earthquakes
Phys. Earth planet. Inter.

200–201
1
9
Frez
J.
Nava
F.A.
Acosta
J.
2010
Source rupture plane determination from directivity Doppler effect for small earthquakes recorded by local networks
Bull. seism. Soc. Am.

100
289
297
Frohlich
C.
1989
The Nature of deep-focus earthquakes
Ann. Rev. Earth Planet. Sci.

17
227
254
Frohlich
C.
2006
Deep Earthquakes

2nd ed.
Cambridge Univ. Press
Fukao
Y.
Obayashi
M.
Inoue
H.
Nenbai
M.
1992
Subducting slabs in the mantle transition zone
J. geophys. Res.

97
4809
4822
Gaherty
J.B.
Lay
T.
Vidale
J.E.
1991
Investigations of deep slab structure using long-period S waves
J. geophys. Res.

96
16 349
16 367
Geman
S.
Geman
D.
1984
Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images
IEEE Trans. Pattern Anal. Mac. Intell.

6
721
741
Glennon
M.A.
Chen
W.-P.
1993
Systematics of deep-focus earthquakes along the Kuril-Kamchatka arc and their implications on mantle dynamics
J. geophys. Res.

98
735
769
N.A.
1964
Total energy and energy spectral density of elastic wave radiation from propagating faults
Bull. seism. Soc. Am.

54
1811
1841
Hayes
G.P.
2011
Rapid source characterization of the 2011 Mw 9.0 off the Pacific coast of Tohoku Earthquake
Earth Planets Space

63
529
534
Hayes
G.P.
et al
2010
Complex rupture during the 12 January 2010 Haiti earthquake
Nature Geoscience

3
800
805
Hayes
G.P.
Wald
D.J.
Johnson
R.L.
2012
Slab1.0: A three-dimensional model of global subduction zone geometries
J. geophys. Res.

117
B01302
doi:10.1029/2011JB008524
Houston
H.
Benz
H.M.
Vidale
J.E.
1998
Time functions of deep earthquakes from broadband and short-period stacks
J. geophys. Res.

103
29 895
29 913
Ingber
L.
1995
Adaptive simulated annealing (ASA): lessons learned
Control Cybern.

25
33
54
Ishii
M.
2011
High-frequency rupture properties of the Mw 9.0 Off-the-Pacific-Coast-of-Tohoku earthquake
Earth Planets Space

63
609
614
Ishii
M.
Shearer
P.M.
Houston
H.
Vidale
J.E.
2005
Extent, duration and speed of the 2004 Sumatra-Andaman earthquake imaged by the Hi-net array
Nature

435
933
936
Kennett
B.L.N.
Engdahl
E.R.
1991
Traveltimes for global earthquake location and phase identification
Geophys. J. Int.

105
429
465
Kennett
B.L.N.
Gorbatov
A.
Spiliopoulos
S.
2014
Tracking earthquake source evolution in 3-D
Geophys. J. Int.

198
867
879
Kikuchi
M.
Kanamori
H.
1991
Inversion of complex body waves–III
Bull. seism. Soc. Am.

81
2335
2350
Kirkpatrick
S.
Gelatt
C.D.
Vecchi
M.P.
1983
Optimization by simulated annealing
Science

220
671
680
Kiser
E.
Ishii
M.
2011
The 2010 Mw 8.8 Chile earthquake: triggering on multiple segments and frequency-dependent rupture behavior
Geophys. Res. Lett.

38
L07301
doi:10.1029/2011GL047140
Kiser
E.
Ishii
M.
2012
Combining seismic arrays to image the high-frequency characteristics of large earthquakes
Geophys. J. Int.

188
1117
1128
Koper
K.D.
Hutko
A.R.
Lay
T.
Ammon
C.J.
Kanamori
H.
2011
Frequency-dependent rupture process of the 2011 Mw 9.0 Tohoku earthquake: comparison of short-period P wave backprojection images and broadband seismic rupture models
Earth Planets Space

63
599
602
Lay
T.
Ammon
C.J.
Kanamori
H.
Rivera
L.
Koper
K.D.
Hutko
A.R.
2010
The 2009 Samoa-Tonga great earthquake triggered doublet
Nature

466
964
968
Lengliné
O.
Got
J.L.
2011
Rupture directivity of microearthquake sequences near Parkfield, California
Geophys. Res. Lett.

38
L08310
doi:10.1029/2011GL047303
Meng
L.
Inbal
A.
Ampuero
J.-P.
2011
A window into the complexity of the dynamic rupture of the 2011 Mw 9 Tohoku-Oki earthquake
Geophys. Res. Lett.

38
L00G07
doi:10.1029/2011GL048118
Meng
L.
Ampuero
J.-P.
Burgmann
R.
2014
The 2013 Okhotsk deep-focus earthquake: rupture beyond the metastable olivine wedge and thermally controlled rise time near the edge of a slab
Geophys. Res. Lett.

41
3779
3785
Metropolis
N.
Rosenbluth
A.W.
Rosenbluth
M.N.
Teller
A.H.
Teller
E.
1953
Equation of state calculations by fast computing machines
J. Chem. Phys.

21
1087
1092
Ni
S.
Kanamori
H.
Helmberger
D.
2005
Energy radiation from the Sumatra earthquake
Nature

434
582
582
Y.
Kasahara
K.
Hori
S.
Obara
K.
Sekiguchi
S.
Fujiwara
H.
Yamamoto
A.
2004
Recent progress of seismic observation networks in Japan: Hi-net, F-net, K-NET and KiK-net
Earth Planets Space

56
xv
xxviii
Pankow
K.L.
Lay
T.
1999
Constraints on the Kurile slab from shear wave residuals sphere analysis
J. geophys. Res.

104
7255
7278
Persh
S.E.
Houston
H.
2004
Deep earthquake rupture histories determined by global stacking of broadband P waveforms
J. geophys. Res.

109
B04311
doi:10.1029/2003JB002762
Pino
N.A.
Mazza
S.
Boschi
E.
1999
Rupture directivity of the major shocks in the 1997 Umbria-Marche (central Italy) sequence from regional broadband waveforms
Geophys. Res. Lett.

26
2101
2104
Pro
C.
Buforn
E.
Udías
A.
2007
Rupture length and velocity for earthquakes in the Mid-Atlantic Ridge from directivity effect in body and surface waves
Tectonophysics

433
65
79
Schwartz
S.Y.
Lay
T.
Beck
S.L.
1991
Shear wave travel time, amplitude, and waveform analysis for earthquakes in the Kurile slab: constraints on deep slab structure and mantle heterogeneity
J. geophys. Res.

96
14 445
14 460
Stein
S.
et al
2011
Learning from failure: the SPREE mid-continental rift experiment
GSA Today

21
5
7
Suzuki
M.
Yagi
Y.
2011
Depth dependence of rupture velocity in deep earthquakes
Geophys. Res. Lett.

38
L05308
doi:10.1029/2011GL046807
Takenaka
S.
H.
Yoshioka
S.
1999
Velocity anomalies and spatial distributions of physical properties in horizontally lying slabs beneath the northwestern Pacific region
Phys. Earth planet. Inter.

112
137
157
Tibi
R.
Wiens
D.A.
Inoue
H.
2003
Remote triggering of deep earthquakes in the 2002 Tonga sequences
Nature

424
921
925
Vallée
M.
Charléty
J.
Ferreira
A.M.G.
Delouis
B.
Vergoz
J.
2011
SCARDEC: a new technique for the rapid determination of seismic moment magnitude, focal mechanism and source time functions for large earthquakes using body wave deconvolution
Geophys. J. Int.

184
338
358
van der Hilst
R.D.
Engdahl
E.R.
Spakman
W.
1993
Tomographic inversion of P and pP data for aspherical mantle structure below the northwest Pacific region
Geophys. J. Int.

115
264
302
Vidale
J.E.
Houston
H.
1993
The depth dependence of earthquake duration and implications for rupture mechanisms
Nature

365
45
47
Wald
D.J.
Kanamori
H.
Helmberger
D.V.
Heaton
T.H.
1993
Source study of the 1906 San Francisco earthquake
Bull. seism. Soc. Am.

83
981
1019
Warren
L.M.
Shearer
P.M.
2006
Systematic determination of earthquake rupture directivity and fault planes from analysis of long-period P-wave spectra
Geophys. J. Int.

164
46
62
Warren
L.M.
Silver
P.G.
2006
Measurement of differential rupture durations as constraints on the source finiteness of deep-focus earthquakes
J. geophys. Res.

111
B06304
doi:10.1029/2005JB004001
Wei
S.
Helmberger
D.
Zhan
Z.
Graves
R.
2013
Rupture complexity of the Mw8.3 Sea of Okhotsk earthquake: rapid triggering of complementary earthquakes?
Geophys. Res. Lett.

40
5034
5039
Wessel
P.
Smith
W.H.
1991
Free software helps map and display data
EOS, Trans. Am. geophys. Un.

72
441
446
Xu
Y.
Koper
K.D.
Sufri
O.
Zhu
L.
Hutko
A.R.
2009
Rupture imaging of the Mw 7.9 12 May 2008 Wenchuan earthquake from back projection of teleseismic P waves
Geochem. Geophys. Geosyst.

10
Q04006
doi:10.1029/2008GC002335
Yagi
Y.
Kikuchi
M.
2000
Source rupture process of the Kocaeli, Turkey, earthquake of August 17, 1999, obtained by joint inversion of near-field data and teleseismic data
Geophys. Res. Lett.

27
1969
1972
Yagi
Y.
Nakao
A.
Kasahara
A.
2012
Smooth and rapid slip near the Japan Trench during the 2011 Tohoku-oki earthquake revealed by a hybrid back-projection method
Earth planet. Sci. Lett.

355–356
94
101
Ye
L.
Lay
T.
Kanamori
H.
Koper
K.D.
2013
Energy release of the 2013 Mw8.3 Sea of Okhotsk earthquake and deep slab stress heterogeneity
Science

341
1380
1384
Zhan
Z.
Kanamori
H.
Tsai
V.C.
Helmberger
D.V.
Wei
S.
2014
Rupture complexity of the 1994 Bolivia and 2013 Sea of Okhotsk deep earthquakes
Earth planet. Sci. Lett.

385
89
96
Zhang
Y.
Feng
W.
Xu
L.
Zhou
C.
Chen
Y.
2009
Spatio-temporal rupture process of the 2008 great Wenchuan earthquake
Science in China Series D: Earth Sciences

52
145
154
Zhou
H.-W.
Anderson
D.L.
Clayton
R.W.
1990
Modeling of residual spheres for subduction zone earthquakes. 1. Apparent slab penetration signatures in the NW Pacific caused by deep diffuse mantle anomalies
J. geophys. Res.

95
6799
6827

## SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this paper:

Table S1. Summary of subevents from back-projection analysis of the 2012 earthquake.

Table S2. Summary of subevents from back-projection analysis of the 2013 earthquake.

Figure S1. Bootstrapping results for Event 1. (a) Histogram of source duration τ obtained from bootstrapping a thousand randomly sampled data sets (grey bars). Thick solid line represents the value obtained from the inversion using all the data. (b) Same as (a) except for the parameter k. (c) Same as (a) except for the rupture dip γ. (d) Same as (a) except for the rupture azimuth ϕ.

Figure S2. Same as Fig. S1 but for Event 2.

Figure S3. Station coverage for the back-projection analyses. (a) Stations (red triangles) used for the analysis of the Mw 7.7 2012 August 14 earthquake. The stations are parts of Alaska Tsunami Warning Seismic System (NOAA/NWS/West Coast-Alaska Tsunami Warning Centre), ANZA Regional Network (IGPP, University of California, San Diego), Berkeley Digital Seismograph Network (Berkeley Seismological Laboratory), Caltech Regional Seismic Network (Caltech/USGS), Canadian National Seismograph Network (Geological Survey of Canada), Global Seismograph Network (IRIS), International Miscellaneous Stations, Lamont-Doherty Cooperative Seismographic Network (Lamont-Doherty Earth Observatory of Columbia University), Leo Brady Network, Ozark Illinois, Indiana and Kentucky ES Flexible Array, United States National Seismic Network (ANSS Data Collection Center), University of Oregon Regional Network (University of Oregon), University of Utah Regional Network (University of Utah), and USArray Transportable Array (IRIS and EarthScope). (b) Same as in (a) except for the Mw 8.3 2013 May 24 earthquake. The blue triangles north of Lake Superior are the additional stations from the SPREE project (Stein et al.2011).

Figure S4. Back-projection result of the 2013 Mw 8.3 earthquake using raw data. Same as in Fig. 11(b) except that the curves shown correspond to back-projection of raw seismograms (i.e. no filtering). Because the data contain lower frequency information, the curves appear smoother than in Fig. 11(b). (http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggv352/-/DC1).

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.