**4. Method introduction**

To perform source parameters studies on the *M*<sup>W</sup> 7.9 mainshock and its larger aftershocks we used several advanced methods. These methods are introduced briefly in this section.

#### **4.1 Method used to obtain a full moment tensor**

An earthquake source can be described using a seismic moment tensor. The moment tensor can be decomposed into three parts: an isotropic (ISO), double-couple (DC), and compensated linear vector dipole (CLVD) part (e.g., [10]). It can also be decomposed into an isotropic (ISO), a major double-couple, and a minor double-couple part (e.g., [11]). The seismic moment is a 3 3 matrices. In linear algebra, a complex matrix can be expressed by the summation of several simple, independent matrices. Applying this principle, Kikuchi and Kanamori [12] expressed an arbitrary moment tensor by summing six different constant moment tensors. Given an earthquake hypocenter and earth model, each of the constant tensors is used to generate Green's functions and obtain three-component synthetic seismograms at a given seismic station.

There are several ways to generate Green's functions depending on the wave type used. As the *M*<sup>W</sup> 7.9 Rat Islands earthquake was large, we used a long-period Rayleigh wave fitting method to obtain its moment tensor. Green's functions were generated using the normal modes summation method [13].

Once Green's functions are obtained, synthetic seismograms are calculated using a set of coefficients. A moment tensor inversion is to search for a set of the coefficients used to generate synthetic seismograms, which should be as similar as possible to the observed seismograms in shapes and amplitudes. To do this, two functions are often separately used. One function is to calculate the correlation between the synthetic and the observed seismograms using:

$$\sigma\_{j} = 1 - \frac{\sum\_{i} X\_{j}(t\_{i}) \bullet O\_{j}(t\_{i})}{\sqrt{\sum\_{i} X\_{j}^{2}(t\_{i}) \bullet \sum\_{i} O\_{j}^{2}(t\_{i})}} \tag{1}$$

where the subscript *j* is an ordinal number of the digital recording; *i* is the data time point index in the observed or the synthetic seismograms. *O*j(*t*i) is a segment of a digital record; *X*j(*t*i) is a segment of a synthetic seismogram corresponding to *O*j(*t*i).

The second function is to calculate the amplitude differences between the synthetic and the observed seismograms:

$$e\_j = \sqrt{\frac{\sum\_i \left[X\_j(t\_i) \times a\_0 - O\_j(t\_i)\right]^2}{N\_j}} \tag{2}$$

The factor *a*<sup>0</sup> is a constant, determined using the following function:

$$\mathfrak{a}\_0 = \frac{1}{N} \sum\_{j=1}^{N} \frac{O\_j(t)\_{\text{max}}}{X\_j(t)\_{\text{max}}} \times 10^{20} \text{ } dyne \text{-} cm \tag{3}$$

where *N* is the number of records used in the inversion.

We, first, use function (1) to obtain a preliminary moment tensor solution, then use (3) to obtain *a*0, and use (2) to obtain a solution. The procedure is repeated many times to find a target set of coefficients when function (2) is at its minimum. Then the moment tensor is calculated using the target set of coefficients. We used the procedures developed by Ma and Adams [14] for simultaneous waveform shape and amplitude inversion.

#### **4.2 The method to determine focal depth using arrival time difference pP-P**

Crotwell, et al. [15] developed a Taup Toolkit, called "Flexible Seismic Travel-Time and Raypath Utilities." Using this tool, the travel times for many seismic phases can be calculated.

We developed a procedure [16] using the relationship that the time duration between the tele-depth phase pP and its reference phase P is roughly positively proportional to the focal depth, to determine focal depths for the larger aftershocks of the Rat Islands, *M*<sup>W</sup> 7.9 earthquake.


### **4.3 Earthquake locating method**

Earthquake hypocenter parameters are fundamental information for studying earthquakes, as such many people have contributed to earthquake locating methods and computer programs. The hypocenter locating program that we used is a part of a computer program package called SEISAN. The SEISAN (seismic analysis system) is a complete set of programs for analyzing earthquakes. With SEISAN it is possible to locate events, determine spectral parameters, seismic moments, and so on. The hypocenter locating program used in this article in the SEISAN package is a modified version of HYPOCENTER [17–19].

#### **4.4 The method to set up a source rupture model**

The commonly used procedure to set up an earthquake rupture model is described below. One of the nodal planes obtained from a seismic moment tensor is used as the earthquake rupture plane. Usually, the *x*-axis is along the strike direction, and the *y*-axis is along the dip direction. The selected rupture plane is divided into *M* � *N* sub-faults with lengths of *dx* and *dy*. Each sub-fault is treated as a point source and the synthetic seismogram at each seismic station is the summation of the synthetic seismograms generated by all of the sub-faults. The source time function of a sub-fault is usually depicted as overlapped triangles. The layout of the source rupture model can be found in Hartzell and Heaton [20].

A unit constant rupture slip vector for each sub-fault is divided into two orthogonal vector components (one aligned along the strike and the other aligned along the dipping direction). Any slip vector on the sub-fault is obtained by multiplying the two constant vector components with appropriate coefficients. The goal of the inversion method is to obtain the coefficients of all of the sub-faults. The slip function (source time function) of each sub-fault is depicted by overlapped *L* triangles with a rise time *τ*, which is half the length of the bottom side of the triangle. The initial constant unit slip direction (slip0) of each sub-fault is the slip of the selected nodal plane. The initial slip is separated into two components in the directions of (slip0 + 45°) and (slip0 - 45°). This breakdown is convenient for Green's functions calculations.

If we assume that on a sub-fault *mn* (*<sup>m</sup>* = 1, ���, *<sup>M</sup>*; *<sup>n</sup>* = 1, ���, *<sup>N</sup>*) at the *<sup>k</sup>*th component direction (*k* = 1, 2), the slip corresponding to the *l* th triangle is *Xmnlk* and the vertical component of the Green's functions generated at station *j* is *gmnkj*, at the same station the vertical component of the synthetic seismogram, *Wj*, generated by all sub-faults at the time point *t*i, is expressed as:

$$\mathcal{W}\_{j}(t\_{i}) = \sum\_{mnlk} X\_{mnlk} \mathbf{g}\_{mnkj} (t\_{i} - (l - 1)\tau - T\_{mn} - dl\_{mn}) \tag{4}$$

where *T*mn is the rupture start time at the *mn*th sub-fault and *dl*mn is a time delay generated by the different travel path lengths between the *P*-waves generated by the *mn*th sub-fault and the rupture start sub-fault *m*0*n*<sup>0</sup> (hypocenter). The set of *Xmnlk* that can generate a *Wj* , which is most similar to the observed seismogram at station *j*, is the best fitting rupture model for the earthquake.

In this study, two methods were used to determine the rupture slip distribution the non-negative least squares (NNLS) method [21] and the simulated annealing (SA) method [22]. For most trial inversions, we used the NNLS method, while the SA method was used at the final step to confirm the solution obtained with NNLS.

The smoothness constraint of the total spatial slip distribution was implemented by a Laplacian differential operator to stabilize the slip solution [23]. To calculate the time delay *dl*mn in function (5), a rupture velocity is required. To calculate Green's functions, we need an initial focal depth, which is also required to obtain a reasonable slip distribution.
