Simultaneous rupture on conjugate faults during the 2018 Anchorage, Alaska, intraslab earthquake (MW 7.1) inverted from strong-motion waveforms

An MW 7.1 ~ 50-km-deep intraslab earthquake within the Pacific/Yakutat slab underlying the North American Plate struck Anchorage, southern Alaska, on November 30, 2018. The ground-motion records very close to the source region of the Anchorage earthquake provide an important opportunity to better understand the source characteristics of intraslab earthquakes in this subduction zone. We estimated the kinematic rupture process during this earthquake using a series of strong-motion waveform (0.05–0.4 Hz) inversions. Our inversions clearly indicate that the Anchorage earthquake was a rare intraslab event with simultaneous rupture on two conjugate faults, which are recognized sometimes for shallow crustal earthquakes but rarely for deep intraslab earthquakes. Interestingly, one of the conjugate faults had low aftershock productivity. This fault extends to great depth and may reflect a deep oceanic Moho or a local low-velocity and high-VP/VS zone within the oceanic mantle. Even though the Anchorage earthquake was a rare event due to the conjugate faults, we found that its kinematic source parameters such as the slip amplitude and large slip area nearly equal the global averages derived from source scaling relationships for intraslab earthquakes. Because the source parameters comparable to the global averages were also found for another large intraslab earthquake in the subducting Pacific/Yakutat slab, these source parameters are likely an important source characteristic common to this subduction zone.


Introduction
At 17:29 Coordinated Universal Time (UTC) on November 30, 2018, an M W 7.1 earthquake occurred in southern Alaska ~ 20 km north of Anchorage at a depth of ~ 50 km (United States Geological Survey (USGS 2019); Fig. 1). Under southern Alaska, the Yakutat terrane, which has characteristics of an oceanic plateau (e.g., Christeson et al. 2010), and the Pacific Plate are moving at the same slip rate but in slightly different directions with respect to the overlying North American Plate (Elliott et al. 2010). The hypocenter of the 2018 Anchorage earthquake was within the subducting Pacific or Yakutat slabs, the exact location of whose interface has not yet been completely defined (e.g., Kim et al. 2014). This hypocenter was also located in the upper part of the transition area from shallow-dipping to steep-dipping subduction (Hayes et al. 2018), and its focal mechanism indicates normal faulting (e.g., Ekström et al. 2012). Therefore, there is no doubt that the Anchorage earthquake was an intraslab event.
Southern Alaska and Anchorage recently experienced two large intraslab earthquakes, the 1999 Kodiak Island (M W 7.0; Hansen and Ratchkovski 2001) and 2016 Iniskin (M W 7.1; Mann and Abers 2020) earthquakes (Fig. 1). Compared to these two events, the 2018 Anchorage earthquake caused stronger ground motions resulting in structural damage to downtown Anchorage because of its closer proximity to Anchorage (e.g., West et al. 2020). During the Anchorage earthquake, the near-source seismograms were recorded by the stations of the United States National Strong-Motion (NP) network and the Alaska Regional (AK) network. These seismograms are very useful for investigating the source rupture of the earthquake. This study estimates the kinematic source rupture process by performing waveform inversions of the near-source strong-motion data at the frequency band of 0.05-0.4 Hz.
Kinematic rupture models of the Anchorage earthquake have already been derived from the inversion of Liu et al. (2019) using geodetic, teleseismic, and strongmotion data and that of He et al. (2020) using geodetic and teleseismic data. Their models indicate that a westward steep-dipping fault rupture is preferred over an eastward shallow-dipping rupture. Even though these studies focused on the rupture during the Anchorage mainshock, the fault geometries in their rupture models insufficiently reproduce the spatially complex aftershock pattern (Fig. 2), as detailed later. Based on this pattern, Ruppert et al. (2020) and West et al. (2020) suggested the possibility of a simultaneous rupture on two conjugate faults during the mainshock. The fault geometry in our inversion considers this possibility and decreases the level of disagreement between the mainshock fault geometry and the aftershock distribution.

Conjugate fault geometry based on the relocated aftershock distribution
We determined the fault planes in our inversion based on the mainshock and the aftershock distributions relocated by the double-difference tomography algorithm . As shown in Fig. 2, the Anchorage mainshock was followed by a vigorous aftershock sequence with ~ 1000 aftershocks (M ≥ 2) within the three days following the mainshock. Because these aftershocks exhibit no clear spatiotemporal migration (Additional file 1: Figure S1), the fault planes assumed for the inversion should cover the entire range of these aftershocks. Most of the aftershocks are concentrated north of the mainshock hypocenter at a depth range of 45-55 km. Meanwhile, the mainshock hypocenter (149.9715° W, 61.3342° N) has a depth of 55.7 km. Even though the relocated hypocenter depth of the mainshock is deeper than the depth of 46.7 km estimated by USGS (2019), the depth of ~ 55 km  Hayes et al. (2018). The dashed bold line demarks the approximate boundary of the Pacific Plate and the Yakutat terrane, and two arrows show their plate motions relative to the overlying North American Plate. The triangles and dashed ellipse show the active volcanoes (Alaska Volcano Observatory (AVO) 2016) and the Cook Inlet Basin, respectively. The inset shows the cross-sectional view along the linear profile X-Y Guo et al. Earth, Planets and Space (2020) 72:176 better matches the arrival times of the observed depth phases in the teleseismic waveforms (Liu et al. 2019). The relocated aftershocks form two evident clusters: a southern cluster as marked by a dashed ellipse and a northern cluster as marked by a dotted ellipse in Fig. 2. The southern cluster dips eastward, and the northern cluster dips westward (Fig. 2b). In general, the aftershocks exhibiting a spatially complex pattern reflect a mainshock rupture on two or more fault planes with different directions, which results both in a difference between the moment tensor and first motion mechanisms and in a low percentage of double-couple component in the moment tensor mechanism, as was the case in the 2008 northern Iwate, Japan, intraslab earthquake (e.g., Suzuki et al. 2009). However, for the Anchorage mainshock, the difference between the moment tensor and first motion mechanisms is very small and the percentage of doublecouple component is significantly high, ~ 90% (USGS 2019). These facts suggest a possibility that the Anchorage mainshock was an event with two conjugate faults corresponding to the two aftershock clusters. The two aftershock clusters, which share a similar focal mechanism , also stimulate us to explore this possibility. Meanwhile, there is an alternative possibility where the mainshock occurred on only one of the conjugate faults and then induced the aftershocks to occur along the other of the conjugate faults. Therefore, as explained later, the inversion of the present study considered several possible cases. Based on the cross-sectional view perpendicular to the strike angle (5°) of the moment tensor solution (Fig. 2b), we found that the southern and northern clusters dip with angles of ~ 30° and ~ 60°, respectively. These dip angles are consistent with those of the two nodal planes in the focal mechanisms for the mainshock and aftershocks . The mainshock hypocenter is likely on the plane of the southern cluster rather than on that of the northern cluster. We also checked another cross-sectional view (Additional file 1: Figure S2), which was perpendicular to the apparent strike angle (340°) of the southern cluster in the map view (Fig. 2a), and confirmed that the dip angles in this direction had large gaps relative to those of the two nodal planes in the focal mechanisms for the aftershocks. Based on the directional superiority of Fig. 2 over Additional file 1: Figure S2, our inversion for the Anchorage mainshock assumed two fault planes, one with a strike of 5° and a dip of 30° and the other with a strike of 185° and a dip of 60° (Fig. 3). These strikes and dips agree well with those of the Global Centroid Moment Tensor (GCMT; Ekström et al. 2012) solution (strike = 7°, 189°; dip = 26°, 64°). As explained later, these planes were divided into three segments (S1,  ) within the three days following the Anchorage mainshock (plotted points). a Map view. The large star marks the epicenter of the mainshock. The two small stars mark the epicenters of the aftershocks used for the calibration of the 1D velocity structure models. The dashed and dotted ellipses denote the southern and northern aftershock clusters, respectively. b Cross-sectional view of the aftershocks along the linear profile A-B. The bold line indicates the slab surface S2, and S3) for our inversion strategy. Segments S1 and S2 represent the deep and shallow portions of the eastward shallow-dipping fault, respectively, and segment S3 represents the westward steep-dipping fault. For simplicity, we assumed an identical fault length of 34 km for these segments. The fault width of segment S1 is 18 km, and segments S2 and S3 have fault widths of 16 km. The hypocenter is located on segment S1. We ruled out the possibility of the hypocenter being located at the intersection of the three segments because the shortest distance from the hypocenter to the intersection (3 km) is greater than the relocation error (< 0.5 km) in Ruppert et al. (2020). The rigidity of the fault planes was assumed to be 56 GPa (Eberhart-Phillips et al. 2006, 2019) without a depth dependence.

Observed data
Our inversion selected nine strong-motion stations ( Fig. 4): seven AK network stations and two NP network stations. The station distribution provides good azimuthal coverage except southwest of the source region. Even though there is a southwestern station  Fig. 1), we excluded this station because strong 3D basin effects (e.g., Mann and Abers 2020; Smith and Tape 2020) complicate the phases of seismic waves, enhance the amplitudes even at extremely low frequency (~ 0.1 Hz), and therefore can adversely affect our inversion results. We confirmed that the nine stations were sufficient to ensure a good spatial resolution of our inversion, even though the station coverage on the western side of the source region was slightly poorer than that on the eastern side. The observed three-component acceleration waveforms were integrated to velocity waveforms, band-pass filtered between 0.05 Hz and 0.4 Hz, and resampled with a sampling frequency of 5 Hz. The upper limit of 0.4 Hz was determined to ensure the accuracy of the Green's functions in our inversion. The waveforms lasting 20-30 s from the manually picked P-wave onsets were used as the inversion data.

Inversion method
We used a multi-time-window linear waveform-inversion method (Hartzell and Heaton 1983), where the rupture process was spatially discretized by dividing the source faults into subfaults (if = 1, nf) and temporally discretized on each subfault by representing the slip histories in two orthogonal slip directions (ir = 1, 2) with multiple time windows (it = 1, nt). A synthetic waveform U syn (t) is modeled as follows: where �(t) is a basis function representing a time window; t , the time interval between consecutive time windows; ξ if , the distance from the hypocenter; V R , the rupture velocity triggering the first time window; and G if ,ir (t) , the Green's function. The model parameters of the waveform inversion were slip amplitudes in the two slip directions for each time window of each subfault ( m if ,it,ir ). All three fault segments in our inversion, S1, S2, and S3, were divided into 2 km × 2 km subfaults. Three time windows, each of which was a smoothed ramp function with a duration of 2 s separated from each other by 1 s were used on each subfault. These time windows permit a maximum slip duration of 4 s and are enough to represent the total rupture on each subfault. These model parameters were obtained by limiting the rake angle (slip direction) to be within ± 45° centered at 270° (pure normal faulting) and solving a non-negative least-squares problem (Lawson and Hanson 1974). For the inversion stability, a spatiotemporal smoothing constraint (e.g., Sekiguchi et al. 2000) was also imposed on the model parameters for each time window: where x and y represent the along-strike and alongdip directions, respectively; and m is the model vector. Therefore, the observation equation is given as where G and d are the matrices comprising the Green's functions and data vector, respectively. The appropriate weight ( ) of the smoothing constraint was selected by referring to the Akaike Bayesian information criterion (ABIC; Akaike 1980; Additional file 1: Figure S3a). The constant V R on each subfault was determined based on (1) the best value of the variance reduction (Additional file 1: Figure S3b), which measures the degree of the fit between the observed and synthetic data: where N, 3, and T are the numbers of stations, waveform components, and time steps, respectively; and U obs (t) represents the observed waveform. We assumed 1D layered velocity structure models and calculated the Green's functions for each subfault using the discrete wavenumber method (Bouchon 1981) and the reflection/transmission coefficient matrix method (Kennett and Kerry 1979). For each station, the 1D layered model was extracted from the 3D gridded velocity structure model of Eberhart-Phillips et al. (2006), (2019). The layered model in the present study was constructed by substituting the physical parameters of the non-grid point with those of the grid point closest to it. Because their 3D model primarily imaged deep structures such as the crust and mantle and therefore did not sufficiently include near-surface low-velocity structures such as sedimentary basins, which amplify and complicate the observed low-frequency (< 1-Hz) waveforms (e.g.  Dutta et al. (2007). We determined the thicknesses of these four layers via trial and error such that the phases and amplitudes in the synthetic waveforms (0.1-0.4 Hz) for two moderate-sized aftershocks (M W 5.0 at 7:57 UTC and M W 4.5 at 12:44 UTC on December 1, 2018; Fig. 2 and Additional file 1: Table S1) satisfactorily reproduced the observations (Additional file 1: Figure S4). Double-couple point sources with smoothed ramp functions and focal information based on Ekström et al. (2012) and Ruppert et al. (2020) were assumed to calculate the synthetic waveforms. Q values based on the formula of Brocher (2008) were used for the calculations. The calibrated 1D velocity structure models ( Fig. 4 and Additional file 1: Table S2)  To verify if a simultaneous rupture on two conjugate faults occurred during the Anchorage mainshock, we performed waveform inversions for three cases (Fig. 3).
(4) variance reduction Case A assumes a rupture on a single, eastward shallow-dipping fault (segments S1 and S2), which is similar to the eastward dipping fault model in the supporting information of Liu et al. (2019). Meanwhile, cases B and C assume a conjugate rupture both on the eastward shallow-dipping and westward steep-dipping faults (segment S3). The difference between the two cases is the absence (case B) or presence (case C) of the shallow portion of the eastward shallow-dipping fault (segment S2). In cases B and C, we assume that the rupture on segment S3 radiated without a time delay from the intersection point (the black star in Fig. 3), with the rupture propagating on segment S1 arriving first. As mentioned earlier, the intersection point should not be the same as the hypocenter, else segment S3 will not agree with the location of the northern aftershock cluster. We also assumed that the slip between the westward steep-dipping (segment S3) and eastward shallow-dipping (segments S1 and S2) faults could be discontinuous. These assumptions follow the inversions of earthquakes with multiple fault planes (e.g., Suzuki et al. 2009). Segments S2 and S3 and the up-dip portion of segment S1 were assumed based on the aftershock distribution. Meanwhile, in all three cases, the down-dip extension of segment S1 (Fig. 3), which resulted in a significantly larger depth relative to the aftershock clusters, was determined via a preliminary inversion trial. Figure 5 shows that, if segment S1 has a small fault width (6 km) that does not extend into the eastern zone with low aftershock activity, an extremely large slip is estimated at its eastern edge. This means that segment S1 should be further expanded. We tested several fault widths of segment S1 (Additional file 1: Figure S5) and determined its width to be 18 km, as mentioned earlier.
Inversion results for the three cases Table 1 summarizes our inversion results, Figs. 6 (cases A, B, and C), 7a (case A), 8a (case B), and 9a (case C) show the final slip distributions, and Figs. 7b (case A), 8b (case B), and 9b (case C) show the distributions of the maximum slip velocities and moment rate functions. The optimal rupture velocities triggering the first time windows for all three cases were determined to be 2.9 km/s (Additional file 1: Figure S3b), which corresponded to 69% of the average S-wave velocity around the source faults. The total seismic moments are 1.2-1.4 times as large as that of the GCMT solution (4.8 × 10 19 Nm). The seismic moment on segment S1 accounts for 61%, 59%, and 44% of the total seismic moment for cases A, B, and C, respectively. The seismic moment of only segment S2 and/or S3 is significantly smaller than that of the GCMT solution. Even though segment S1 had low aftershock activity that led us to debate whether the mainshock rupture extended to it, the large percentage of segment S1   suggests that it ruptured during the mainshock. Furthermore, if segment S1 had small fault width shown in Fig. 5, the observed waveforms at the eastern stations away from the source region could not be reasonably reproduced (Additional file 1: Figure S5b). The spatial distributions of the estimated large slips and aftershocks are complementary to each other (Fig. 6). On segments S1 and S2 (cases A and C), large slips occur east (~ 60 km in depth) and north (~ 50 km in depth) of the hypocenter, even though the northern slip in case C is slightly uncertain. These large slips surround most of the southern aftershock cluster. On segment S3 (cases B and C), the large slip is located at a depth of ~ 50 km just south of the densest clustering of northern aftershocks. We can directly compare the slip distribution of case A with that of the eastward shallow-dipping model in Liu et al. (2019) because the two have similar fault geometries. The slip distribution of case A is consistent with that of Liu et al. (2019) in that the large slips are located a small distance from the hypocenter.
In case B, the peak values of the large slips on the conjugate segments S1 and S3 are comparable (Figs. 6 and 8a). In case C, the added segment S2 results in a slip separation from S3 to S2 and therefore in a reduction of the large slip on S3. The spatial slip patterns on segments S2 and S3 in case C are similar to those in cases A and B, respectively. In case C, even though the large slip on segment S3 is reduced, its peak value is still ~ 1.5 times greater than the maximum slip on segment S2 (Figs. 6 and 9a).

Determination of the preferable case
Cases A and B have comparable total variance reductions (Table 1). However, we favor case B, that is, the westward steep-dipping fault in the upper half of the source (segment S3), based on its better reproduction of the observed main pulses at station AHOU, which lies adjacent to the northern edge of the source fault and recorded considerable amplitudes. Two large pulses can be seen at a lapse time of 10-15 s in the east-west component (Fig. 10b). Segment S3 in case B reproduces the amplitudes of both pulses reasonably well, even though segment S2 in case A reproduces only the second pulse. This is primarily due to the upward rupture propagation on the steep-dipping segment S3. To confirm that case B was better than case A, we used two virtual line source models and performed a numerical test (Fig. 11). The rupture on line source models A and B propagated toward the farthest point on segments S2 and S3, respectively. Note that in this numerical test, segment S1 was not considered because it is a segment common to inversion cases A and B. For calculating the synthetic waveforms, we assumed a unit slip, pure normal faulting, and a rupture velocity of 2.9 km/s. The spatial pattern of the synthetic peak amplitudes for line source A, where the peak amplitude for the station on the western side of the hypocenter (SSN) was larger than that for the stations on the eastern side and above the thick sediments (AHOU and AWCH; Fig. 4), cannot be found in the observation. This fact suggests that it is likely that the main rupture did not extend into segment S2 in inversion case A. On the other hand, the synthetic spatial pattern for line source B, which includes the northeastward rupture propagation on the steep-dipping segment S3, is similar to the observation. Analyses of interferometric synthetic aperture radar images (e.g., He et al. 2020;West et al. 2020) also indicate that the westward steep-dipping fault is required to explain the significant uplift on the eastern side of the source region.
Together with the comparison between cases A and B (Figs. 10b and 11), the comparison between cases B and C also eliminates the possibility of the mainshock rupture on segment S2. Table 1 shows that the variance reduction of case C is significantly higher than that of case B, which means that the total waveform fit (Fig. 10) is improved by the addition of segment S2. However, this improvement is not surprising because the added segment S2 increases the number of model parameters and therefore allows a more detailed waveform fit. In general, a model with more degrees of freedom and a higher variance reduction is not always better in that it may not correspond to other datasets. To judge whether case B or C is better, we calculated the values of the Akaike information criterion (AIC; Akaike 1974), which considers both the data fit and the model complexity/simplicity and is often used for model selection. We found that case B (2.0 × 10 2 ) has a smaller AIC value than case C (9.7 × 10 2 ). Therefore, we suggest that case B is the better model and that most of the actual rupture during the Anchorage mainshock was restricted to segments S1 and S3 and did not extend to segment S2. Even though the complex aftershock distribution and the large focal depth made it difficult to precisely determine the extent of the source faults during the Anchorage mainshock, the above comparisons lead us to conclude that case B is the most plausible of the three cases.

Discussion on the source characteristics of case B
Here, we detail the source characteristics of case B, our best model. Because large slips have a major influence  Somerville et al. (1999) on ground motions, an understanding of the characteristics of the large slips during the Anchorage earthquake will contribute to ground-motion predictions of future intraslab earthquakes. To identify the large slip area in case B, we used the criteria proposed by Somerville et al. (1999) as follows. First, we decided that the rupture area for case B was the same as the combined area of the modeled segments S1 and S3 because the average slips on the fault-edge rows/columns for case B (Fig. 8a) were greater than 0.3 times the average slips over the modeled segments. Then, we defined the rectangular area of the subfaults with ≥ 1.5 times the average slip of the rupture area as a large slip area (called "asperity area" in Somerville et al. (1999)). As a result, we obtained two large slip areas in case B. One (A1) is on the down-dip side of the hypocenter on segment S1, and the other (A2) is in the central portion of segment S3 (Fig. 8). A1 and A2 were generated 2.5-5 s and 5-10 s after the rupture initiation, respectively (Fig. 12a). As illustrated in the contributions of A1 and A2 to the velocity waveforms at three near-source stations (Fig. 12b), much of the amplitude can be explained by these large slip areas and the contribution of A2 occurs subsequent to that of A1. The synthetic peak amplitudes are produced by A2, and its value at the northern station AHOU is the largest due to the forward directivity resulting from the upward and northward rupture propagation of A2. At stations AHOU and K209, A2 yields larger amplitudes than A1 because of the forward directivity and the shallower depth of A2, respectively. Meanwhile, the amplitude from A1 is comparable to that from A2 at the eastern station K217 because of the eastward rupture propagation of A1. The two large slip areas have a total dimension that is 17% of the rupture area. This percentage is equivalent to that calculated by Iwata and Asano (2011), who proposed scaling relationships of the source parameters for M 7-8-class intraslab earthquakes. Even though several global source scaling relationships for intraslab earthquakes have been published (e.g., Iwata and Asano 2011;Allen and Hayes 2017), these relationships were developed using no source parameter data from southern Alaskan (i.e., Alaska-Aleutian) earthquakes because of the low intraslab activity in this area compared to other subduction zones. In comparison with Iwata and Asano (2011), we found that the large slip areas of the Anchorage earthquake had a globally average total dimension (Fig. 13a).
The total dimension and average slip of the rupture area are also key parameters for constructing source models of scenario earthquakes. We compared these parameters of the Anchorage earthquake to source scaling relationships (Iwata and Asano 2011;Allen and Hayes 2017). Such scaling relationships can also account for the dimension (Fig. 13b) and average slip (Fig. 13c) of the rupture area of the Anchorage earthquake, as well as the dimension of the large slip area. Figure 13 also compares the Anchorage earthquake to the 2016 Iniskin earthquake (USGS 2018) in southern Alaska, the source parameters (Additional file 1: Figure S6) of which were also obtained using the criteria of Somerville et al. (1999). The previously described feature is shared by the Iniskin earthquake. Even though source parameters include regional perturbations associated with differences in the subduction environment (e.g., Stirling et al. 2013), we suggest that the source parameters of intraslab earthquakes in southern Alaska nearly equal the global average source parameters.
During the Anchorage earthquake, the large slip area A2 propagated in the direction away from downtown Anchorage (Fig. 12). Despite the source parameters comparable to the global averages and the backward direction, ground motions were observed with large low-frequency amplitudes (e.g., Moschetti et al. 2020).

Fig. 10
Comparison between the observed and synthetic velocity waveforms (0.05-0.4 Hz) during the mainshock. a For each station, the observed waveforms and synthetic waveforms for cases A, B, and C are arranged in order from top to bottom. The amplitudes of each waveform were normalized to the peak amplitudes of the observed three-component waveforms. The number to the right of each waveform denotes the peak amplitude in cm/s. b Detailed view of the east-west component of station AHOU. The two arrows indicate the main pulses. The observed waveform, synthetic waveforms calculated by assuming segments S1 and S2 (case A) and segments S1 and S3 (case B), and synthetic waveforms calculated by assuming only segment S2 (case A) and segment S3 (case B) are arranged in order from top to bottom Figure 14 shows the velocity Fourier spectra at several stations in downtown Anchorage. The eastern stations (ALUK and K215), which are located near the edge of the Cook Inlet Basin, have no clear spectral peaks in the low frequency. This means that the observed large low-frequency ground-motions did not originate from the earthquake source. Meanwhile, the western stations farther from the basin edge (8036 and K220) have remarkably large amplitude levels at < 0.5 Hz and spectral peaks at 0.2-0.3 Hz. The spectral difference between the eastern and western stations suggests that the large low-frequency ground-motions in downtown Anchorage were primarily due to site amplification effects associated with the basin, rather than source effects.

Tectonic structure around the source region
The oceanic crust is generally characterized by a lowvelocity and high-V P /V S zone. Because this zone suggests the existence of the fluid released by dehydration reactions, which reduces the effective normal stress and fault strength, intraslab earthquakes usually occur within this zone (e.g., Raleigh and Paterson 1965). Around the source region of the Anchorage earthquake, a laterally extended zone with low-velocity and high-V P /V S (> 1.8) was detected at a depth of ~ 50 km (e.g., Eberhart-Phillips et al. 2006Kim et al. 2014; Additional file 1: Table S2). This means that the oceanic Moho of the subducting Pacific/Yakutat slab extends to at least ~ 50 km. As shown in Fig. 3, the lower limit of our conjugate fault model penetrated to a greater depth of 63 km. Note that, due to the eastward shallow-dipping fault (segment S1), this depth was shallower than the ~ 70 km found by Liu et al. (2019) and He et al. (2020), whose best models were on a single, westward steep-dipping fault. Our conjugate fault model suggests two possibilities. The first is that the oceanic Moho has a depth of 60-65 km at the edge of the shallow-dipping Pacific/Yakutat slab. If this is true, the total thickness of the oceanic crust from the slab surface (30-35 km) is ~ 30 km. The second possibility is that, even though the oceanic Moho is located at ~ 50 km, there is a local low-velocity and high-V P /V S zone within the oceanic mantle (e.g., Nakajima et al. 2011) around the source region of the Anchorage earthquake. To distinguish between these possibilities, it is important to conduct further investigations of seismic structure imaging with finer resolution. For each station, the first waveforms from the top show a comparison between the observation and the synthetic calculated for the entire fault. The second and third waveforms from the top show the synthetics calculated using only the large slip areas A1 and A2, respectively. The map view shows the location of the stations (triangles), the epicenter (star), and the large slip areas (rectangles)

Conclusions
We performed strong-motion waveform (0.05-0.4 Hz) inversions and proposed a kinematic rupture model of the 2018 Anchorage intraslab earthquake (M W 7.1) based on the aftershock distribution as relocated by Ruppert et al. (2020). Our inversions revealed a simultaneous rupture on two conjugate faults, an eastward shallowdipping fault and a westward steep-dipping fault. The preferred rupture model suggests that the shallower side of the eastward shallow-dipping fault, where aftershocks formed a distinct cluster, experienced little rupture during the mainshock. Our conjugate faults also showed the possibility of a deep oceanic Moho or a local low-velocity and high-V P /V S zone within the oceanic mantle of the subducting Pacific/Yakutat slab. On each of the conjugate faults, we identified one large slip area that significantly contributed to the observed ground motions. Unlike shallow crustal earthquakes, there have been few identified large intraslab earthquakes with conjugate faults. We found that, even though the Anchorage earthquake was a rare intraslab event in terms of having conjugate faults, its source parameters, such as its slip amplitude and fault dimension, are consistent with the averages of global intraslab earthquakes.
Additional file 1: Table S1. Focal Mechanism of Two Moderate-Sized Aftershocks. Table S2. Physical Parameters in the 1D Velocity Structure Models Used for the Waveform Inversions. Figure S1. Spatiotemporal distribution of the relocated aftershocks within the three days following the Anchorage mainshock (plotted points). The large star marks the epicenter of the mainshock. Figure S2. Relocated aftershocks within the three days following the Anchorage mainshock (plotted points). a Map view. The large star marks the epicenter of the mainshock. The dashed and dotted ellipses denote the southern and northern aftershock clusters,