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.
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.
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)
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.
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
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 = L/τV, 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,
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,
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
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,
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,
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).
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).
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).
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
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
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.
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.
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).
|Event||τ, t (s)||k||γ (deg)||ϕ (deg)||VR (km s−1)||L (km)|
|1||26 ± 1||0.27 ± 0.02||48 ± 8||42 ± 5||2.7 ± 0.2 (49 ± 4 per cent)||69 ± 6|
|2||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−1)||L (km)|
|1||26 ± 1||0.27 ± 0.02||48 ± 8||42 ± 5||2.7 ± 0.2 (49 ± 4 per cent)||69 ± 6|
|2||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.
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).
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.
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 (tA − tB) 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.
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,
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.
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.
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.