Consecutive ruptures on a complex conjugate fault system during the 2018 Gulf of Alaska earthquake

We developed a flexible finite-fault inversion method for teleseismic P waveforms to obtain a detailed rupture process of a complex multiple-fault earthquake. We estimate the distribution of potency-rate density tensors on an assumed model plane to clarify rupture evolution processes, including variations of fault geometry. We applied our method to the 23 January 2018 Gulf of Alaska earthquake by representing slip on a projected horizontal model plane at a depth of 33.6 km to fit the distribution of aftershocks occurring within one week of the mainshock. The obtained source model, which successfully explained the complex teleseismic P waveforms, shows that the 2018 earthquake ruptured a conjugate system of N-S and E-W faults. The spatiotemporal rupture evolution indicates irregular rupture behavior involving a multiple-shock sequence, which is likely associated with discontinuities in the fault geometry that originated from E-W sea-floor fracture zones and N-S plate-bending faults.


Introduction
This Supporting Information contains numerical tests for validation of the developed finite-fault inversion method (Material S1, Figs. S1 to S4, and Table S1). Sensitivity of the finite-fault inversion to assumptions of model planar depth and rupture velocity is shown in Material S2, and Figs. S5 and S6. Comparison with the conventional smoothness constraints is shown in Figs. S7 and S8. Results of a numerical test using our solution of the 2018 Gulf of Alaska earthquake as input data are shown in Fig. S9. The possibility of dummy imaging of reverberations is evaluated in Material S3 and Fig. S10. Waveform fits and full snapshots of the rupture evolution for our main result are shown in Figs. S11 and S12. Fig. S13 shows a comparison between our source model and the back-projection result obtained by Krabbenhoeft et al. (2018). Table S2 shows the near-filed velocity structure used for calculating Green's function. Table S3 provides the set of smoothness constraints adopted for the main result.

Numerical experiment (S1)
We perform the numerical tests to evaluate effects of the improved smoothness constraints and the horizontal non-rectangular model plane. To generate synthetic waveforms, orthogonal three faults were assumed (Fig. S1a). Then, we assume the pure strike-slip rupture which spherically spread from a hypocenter at a depth of 30 km on the central fault, named F2, with a rupture velocity of 3.0 km/s ( Fig. S1a and b). The moment rate function of input model has peaks at 9 and 23 s (Fig. S1c). We add a random Gaussian noise to the calculated Green's function, for which the standard deviation is 5% of maximum amplitude of each calculated Green's function. We also add a random Gaussian noise with zero mean and a standard deviation of 1.0 as background noise. We generate the synthetic waveforms at 78 stations used in the inversion of the 2018 Gulf of Alaska earthquake (Fig. 2c).
We compare results of four cases: (1) the rectangular model plane and the conventional smoothness constraints; (2) the rectangular model plane and the improved smoothness constraints; (3) the non-rectangular model plane and the conventional smoothness constraints, and (4) the non-rectangular model plane and the improved smoothness constraints.
We set the horizontal rectangular model plane with a width and length of 120 km to cover the input three faults (Fig. S2a). The depth of the model plane is set to 30 km, which corresponds to the centroid depth of input source model. The spatial knot interval is set to 10 km. For the cases (3) and (4), we design the non-rectangular model plane based on the input three faults (Fig. S2d). For all cases, the potency-rate density function at each knot is represented as a linear combination of B-spline functions over a duration of 30 s with an interval of 0.8 s and the rupture front velocity set at 7.0 km/s. We adopt the improved smoothness constraints at the cases (2) and (4) by referring to the input focal mechanism (Table S1).
In the case (1), the resultant moment rate function is smoother than the input one and has only one peak at 12 s, which is about 3 s later than the first peak of the input (Fig.  S2b). The normalized L2 norm, which represents the degree of misfit between the input and the resultant moment rate function (hereinafter called "the L2 norm" for simplicity), was 0.245. The snapshots show a wider potency-rate density distribution than the input, making it difficult to identify the fault geometry and interpret the source process ( Fig. S3a and b). Figure S4 shows the self-normalized potency-rate function for each basis component, obtained by taking a spatial integration of the potency-rate density function for each basis component. In the case (1), the potency-rate function of 1 component (Kikuchi & Kanamori, 1991), corresponding to the input slip direction, is smoother than those of other components (Fig. S4). This is because the conventional smoothness constraints work to excessively smooth out the dominant basis components.
In the case (2), the moment rate function yields two peaks at 10 and 23 s, which close to the input peaks (Fig. S2c). The L2 norm is 0.074. The improved smoothness constraints remove the bias in the resultant potency-rate function of the 1 component ( Fig. S4) and thus the spatiotemporal potency-rate density distribution of the case (2) is finer than that of the case (1) (Fig. S3c). However, in the case (2), the image looks too blurry to resolve two independent ruptures of the input model due to insufficient spatial resolution (snapshot at 15 and 20 s in Fig. S3c).
In the case (3), the moment rate function has two peaks at 11 and 23 s, which close to the input peaks (Fig. S2e). The L2 norm is 0.135 and slightly larger than the case (2). The potency-rate density distribution of the case (3) at 15 and 20 s resolves two ruptures, which are not well resolved in the case (2) (15 and 20 s in Fig. S3d). The spatial resolution of the inversion results is improved because the model space modification according to the input fault geometry is identical to implicitly introducing a priori constraint of the fault geometry (e.g., aftershock distribution). The model space reduction also contributes to reduce computational costs, which is useful for analyses of earthquakes having a vast source area, such as the 2018 Gulf of Alaska earthquake.
In the case (4), the moment rate function reproduces the input in detail (Fig. S2f), and the rupture evolution is fine enough to reproduce the input (Fig. S3e). The L2 norm of the case (4) is 0.071, which is the minimum value among the four cases. Thus, we conclude from our numerical tests that the optimum strategy should be by using both the improved smoothness constraints and the horizontal non-rectangular model plane.

Sensitivity test (S2)
We evaluate the sensitivity of the inversion results by perturbating the model parameters. We perform the inversion analyses by changing the model plane depth to 33.6±5 km. The obtained snapshots show the rupture pattern is insensitive to the model planar depths (Fig. S5). We also check the inversion results by changing the assumption of maximum rupture velocity to 3 and 5 km/s. We resolve the similar rupture processes for the maximum rupture velocities at 5 and 7 km/s ( Fig. S6b and c). However, when assuming 3 km/s, the model does not clearly show the A2 rupture (Fig. S6a). This is due to the limited model space that could artificially vanish the possible slip behavior beyond the designated rupture front.

Verification of using empirical Green's functions (S3)
As shown in Figs. 2 and S11, our finite-fault model sufficiently reproduces the complicated observed teleseismic P waves, resulting in showing the complex-multiple rupture episodes. On the other hand, pulses of the observed waveforms may include later arrivals due to structure complexities in the source region (e.g., Fan & Shearer, 2018;Yue et al., 2017). However, the theoretical Green's functions, assuming a 1D-layered structure model, are often poorly modeled for reverberations of dipping near-source bathymetry (Wiens 1987(Wiens , 1989, which may induce artificial imaging of multiple-shock sequence. In principle, seismograms of relatively smaller earthquakes with a similar focal mechanism that occurred near the target earthquake can be regarded as an empirical Green's function (EGF) under an assumption that the moment rate function of that small earthquake is simple and short (Hartzell 1978;Dreger 1994). We here employ the EGFs instead of the theoretical Green's functions to evaluate whether multiple-shock sequence that we resolve is likely from the source effect or the reverberations. We deconvolve the EGFs from the observed waveforms of the 2018 Gulf of Alaska earthquake for each station to remove the effects of the earth response including possible reverberations.
We select three events from the GCMT catalog (Dziewonski et al., 1981;Ekström et al., 2012) as the EGFs with the clear first P-phase motion and high signal-to-noise ratios (Fig. S10a). The EGFs and the mainshock data are band-passed between 0.01 and 2 Hz and converted into ground velocities with a sampling interval of 0.1 s. We solve the least squares problem using the non-negative least squares algorithm of Lawson and Hanson (1974). We perform deconvolution for both a maximum source duration of 65 and 27 s to evaluate the validity of the sub-events resolved after 27 s for the mainshock (Fig. S10b, c,  d, and e).
The normalized moment rate functions obtained in the maximum length of 65 s show non-negligible moment release even after 27 s (Fig. S10b, c, d and e). If the subevents after 27 s were artifacts caused by the reverberations of the initial rupture, the observed waveforms would be reproducible by convolving the moment rate function up to 27 s with the EGF. However, the synthesized waveforms obtained from the 27-s-moment-rate function fails to reproduce the several pulses of the observed waveforms, while the synthesized waveforms obtained from the 65-s-moment-rate function better fits the observed waveforms, suggesting that the subevents after 27 s should be necessary to explain the observed data (Fig. S10b, c, d and e).      (Matthews et al, 2011;Wessel et al., 2015), respectively. The large beachball in each panel indicates the corresponding total moment tensor at each time. This figure was made with matplotlib (v3.1.1; Hunter, 2007) and ObsPy (v1.1.0;Beyreuther et al., 2010).   (Matthews et al, 2011;Wessel et al., 2015), respectively. The large beachball in each panel indicates the corresponding total moment tensor at each time. This figure was made with matplotlib (v3.1.1; Hunter, 2007) and ObsPy (v1.1.0;Beyreuther et al., 2010). Figure S8. Comparison of the potency-rate functions of the 2018 Gulf of Alaska earthquake obtained by the conventional (left column) and improved smoothness constraints (right column) for each basis double-couple component (Kikuchi & Kanamori, 1991). Each trace is self-normalized. This figure was made with matplotlib (v3.1.1; Hunter, 2007) and ObsPy (v1.1.0; Beyreuther et al., 2010).  (Bird, 2003) and the fracture zones (Matthews et al, 2011;Wessel et al., 2015), respectively. This figure was made with matplotlib (v3.1.1; Hunter, 2007) and ObsPy (v1.1.0;Beyreuther et al., 2010). Figure S10. Summary of the EGF analysis. (a) Map projection of the GCMT solutions of the main shock (orange beachball) and events used as the EGFs (black beachballs). The star is the mainshock epicenter, and orange dots are aftershocks (M ≥ 3) that occurred within one week of the mainshock; all epicentral locations are from AEC. Dashed and solid lines represent the plate boundaries (Bird, 2003) and the fracture zones (Matthews et al, 2011;Wessel et al., 2015), respectively. Inset is an azimuthal equidistant projection of the station distribution. (b) to (e) show the normalized moment rate function (left) and waveform fittings (right). Gray trace is the observed waveform. Also shown are the synthetic waveforms obtained by using 65-s-moment-rate function (orange) and 27-s-moment-rate function (blue). Each panel is labeled with the station name, azimuth (Azi.) and epicentral distance (Del.) from the mainshock, and the event name used as the EGF. This figure was made with matplotlib (v3.1.1; Hunter, 2007) and ObsPy (v1.1.0;Beyreuther et al., 2010).  Figure S12. Snapshots of the potency-rate density tensors every 1 s for the 2018 Gulf of Alaska earthquake. The dotted line shows the border of the assumed model plane. The star and solid lines indicate the epicenter (AEC) and the fracture zones (Matthews et al, 2011;Wessel et al., 2015), respectively. The large beachball in each panel indicates the corresponding total moment tensor at each time. This figure was made with matplotlib (v3.1.1; Hunter, 2007) and ObsPy (v1.1.0; Beyreuther et al., 2010). Figure S13. Comparison between the potency-rate density obtained by this study and the high-frequency (0.5-2.0 Hz) emissions obtained by Krabbenhoeft et al. (2018), projected along (a,b) north-south and (c,d) east-west directions from the epicenter. Positive distance directs toward (a,b) north and (c,d) east. Both the potency-rate density and semblance peaks distributions show general agreement in rupture pattern, especially during 0 to 10 s from the origin time. We note that a direct comparison or superimposing is not made here because our study and Krabbenhoeft et al. (2018) adopted the different referencerupture points (epicenters). The dashed line represents the reference rupture speed. This figure was made with matplotlib (v3.1.1; Hunter, 2007). Table S1. Factors of the smoothness constraint of each potency component for the numerical tests. The number represents to components defined by Kikuchi and Kanamori (1991). * * is the absolute value of the total potency derived from the input total moment tensor (Fig. S1b). The scaling factor was set so that ( * *) = .
Table S2. CRUST1.0 structural velocity model (Laske et. al., 2013). Table S3. Factors of the smoothness constraint of each potency component for the analysis of the 2018 Gulf of Alaska earthquake. * * is the absolute value of the total potency of each potency component derived from the GCMT solution (Fig. 1). The scaling factor was set so that ( * *) = .