Focal Mechanisms and the Stress Field in the Aftershock Area of the 2018 Hokkaido Eastern Iburi Earthquake (MJMA = 6.7)

The tectonic stress eld was investigated in and around the aftershock area of the Hokkaido Eastern Iburi earthquake (M JMA = 6.7) occurred on 6 September 2018. We deployed 26 temporary seismic stations in the aftershock area for approximately 2 months and located 1785 aftershocks precisely. Among these aftershocks 818 focal mechanism solutions were determined using the rst motion polarity of P wave from the temporary observation and the permanent seismic networks of Hokkaido University, Japan Meteorological Agency (JMA), and High Sensitivity Seismograph Network Japan (Hi-net). We found that (1) the reverse faulting and the strike-slip faulting are dominant in the aftershock area, (2) the average azimuths of P- and T-axes are N78° ± 33°E and N3° ± 52°W, respectively, and (3) the average dips of P-and T-axes are 25° ± 16° and 46° ± 20°, respectively: the P-axis is close to be horizontal and the T-axis is close to be vertical. We applied a stress inversion method to the focal mechanism solutions to estimate a stress eld in the aftershock area. As a result, we found that the reverse fault type stress eld is dominant in the aftershock area. An axis of the maximum principal stress (σ 1 ) has the azimuth of N73° ± 8°E and the dipping eastward of 17° ± 6° and an axis of the medium principal stress (σ 2 ) has the azimuth of N126° ± 91°E and the dipping southward of 16° ± 13°, indicating that both of σ 1 - and σ 2 -axes are close to be horizontal. An axis of the minimum principal stress (σ 3 ) has the dipping westward of 64° ± 9° that is close to be vertical. The results strongly suggest that the reverse-fault-type stress eld is predominant as an average over the aftershock area which is in the western boundary of the Hidaka Collision Zone. Although the average of the stress ratio is R = 0.6 ± 0.2 in the whole aftershock area, R decreases systematically as the depth is getting deep, which is modeled by a quadratic polynomial of depth.


Introduction
The tectonic regime is complicated in Hokkaido corner, Japan subduction zone (Fig. 1). The Paci c (PA) plate is moving toward N63°W with a speed of 8.2 cm/y (DeMets et al 1994) and subducting below the North American (NA) plate or the Okhotsk (OK) plate on which the Hokkaido Island is located (Takahashi et al 1999;Katsumata et al 2002Katsumata et al ,2003. The upper surface of the PA plate strongly coupled with the overriding plate in and around the Kurile Trench (Hashimoto et al 2009) and shallow great earthquakes have been caused repeatedly. This subduction process possibly produces a compressional stress eld with the direction of NW-SE in the inland area of Hokkaido Island. Moreover, a collision process is in progress. The Kurile Islands arc and the NE Japan arc are colliding in and around the Hidaka Mountain Range (HMR) (Kimura 1981(Kimura ,1986(Kimura ,1996Seno 1985;Moriya 1986;Arita et al 2001). This is called the Hidaka Collision Zone (HCZ). The speed of the collision is estimated to be 6-11 mm/y relative to the NA plate based on the horizontal slip direction from shallow-thrust earthquakes (DeMets 1992). This collision process possibly produces a compressional stress eld in the HCZ with the direction of NE-SW.
Additionally, the upper crust and the lower crust beneath the HCZ are not a simple layered structure (Ozel et  Some tectonic models have been proposed for the HCZ. The crust of the Kurile Islands arc has been torn in the east of the HMR due to the collision and divided into the upper part and the lower part (Ito et al 1999;Murai et al 2003). The upper part is riding over the NE Japan arc in the west of the HMR. The lower part is in contact with the upper boundary of the PA plate, dragged into the upper mantle and scraped (Moriya 1999;Tsumura et al 1999). Kita et al (2012) insisted that the mantle material might be rising directly from the uppermost mantle of the Kurile Islands arc. These complicated structures may cause the complicated stress eld and produce earthquakes with a various type of the focal mechanisms. In the western boundary of the HCZ, a large earthquake occurred on 6 September 2018: The Hokkaido Eastern Iburi earthquake (M JMA = 6.7). Although the focal mechanism of the main shock was estimated as a strike-slip faulting by using the rst motion polarities of P wave (JMA 2018b; NIED 2018b; Katsumata et al 2019), the centroid moment tensor (CMT) solution shows the reverse faulting with the Paxis in the direction of NE-SW (JMA 2018a; NIED 2018a). The mismatch between the focal mechanism solution and the CMT solution has been explained by a model that a large reverse faulting occurred immediately after an initial rupture of a small strike-slip faulting (Katsumata et al 2019). The CMT solution of the 2018 Hokkaido Eastern Iburi earthquake is similar to the focal mechanisms of the 1970 Hidaka earthquake and the 1982 Urakawa-oki earthquake, clearly indicating that the compressional stress eld due to the collision extends to the western boundary of the HCZ.
The main shock of the 2018 Hokkaido Eastern Iburi earthquake was followed by lots of aftershocks. We deployed temporary seismic stations densely in the aftershock area to determine the hypocenters and the focal mechanisms accurately. The purpose of this study is to determine the focal mechanisms of the aftershocks, to apply a stress inversion method to the focal mechanisms, and to make some discussions on the detailed spatial pattern of the stress eld in the aftershock area.

Data
To obtain aftershocks data in detail, we deployed 26 temporary seismic stations in the focal area immediately after the main shock and observed aftershocks for approximately 2 months ( Fig. 2 and Additional le 1). The temporary observation was conducted by the Group for the Aftershock Observations of the 2018 Hokkaido Eastern Iburi Earthquake, which consists of Hokkaido University, Hirosaki University, Tohoku University, Chiba University, the University of Tokyo, Nagoya University, Kyoto University, Kyushu University, Kagoshima University, and the National Research Institute for Earth Science and Disaster Resilience (NIED). The temporary seismographic stations consisted of 4 telemetry online systems and 22 portable o ine systems. We also used 183 permanent online seismographic stations maintained by Hokkaido University, JMA, and NIED. Waveform data observed during the period from 6 September 2018 to 31 October 2018 were examined carefully by visual inspection, and the arrival times of P and S waves and the rst motion polarities of P wave were read manually by a well-trained person.

Methods
We determined hypocenters of earthquakes with the maximum likelihood estimation algorithm of Hirata and Matsu'ura (1987) using the 1D velocity structure of P wave based on Kasahara et al (1994) (Fig. 2), which is the same as that used for the hypocenter calculation at the Hokkaido University. The S wave velocity was obtained by the relationship , where V p and V S are the P and S wave velocities, respectively. We located 1785 earthquakes in the study area ( We determined focal mechanism solutions of earthquakes by using a grid-search technique developed by Hardebeck and Shearer (2002). We used two 1D velocity structures to take the uncertainty of ray paths, especially take-off angle from the hypocenter, into account (Fig. 2). The rst one is a hybrid of two previous studies: the P wave velocity in the crust shallower than 10 km is based on a refraction experiment (Iwasaki et al 2004)  There are several stress inversion methods to estimate the state of stress from focal mechanisms. Only four independent components are able to be obtained from the stress inversion method: the orientation of the axes of three principal stresses and the stress ratio. The principal stresses are the maximum principal stress (σ 1 ), the intermediate principal stress (σ 2 ), and the minimum principal stress (σ 3 ) and the stress ratio R is de ned as R = (σ 1 -σ 2 ) / (σ 1 -σ 3 ), indicating the relative magnitude of the principal stresses and ranging from 0 to 1. We determined the stress eld with a stress inversion method developed by Hardebeck and Michael (2006) using focal mechanisms as input data. The method is performed by dividing the study area su ciently ne in advance, putting a constraint that the stress changes smoothly between neighboring areas to avoid instability of the solution, and calculating the stress of all areas at once by using a least squares method. The uncertainty of the parameters is estimated using 2000 bootstrap resampling of all data (Hardebeck and Michael 2006). In this study the two-dimensional nodes are placed in the aftershock area: the latitude ranges from 42.55 to 42.85°N and the grid spacing is 0.05°, the longitude is xed at 142.0°E for all nodes, and the depth ranges from 8.2 to 45.1 km and the grid spacing is 4.1 km. The focal mechanisms of aftershocks that occurred within 7 km from each node were used. We selected the damping parameter e (Eq. (14) in Hardebeck and Michael 2006) based on the trade-off curve between the model length and the data variance. The corner of the trade-off curve was near e ≈ 1.2, so we selected e = 1.2 for all groups in this study.

Focal mechanisms
We determined 818 focal mechanisms from 1785 aftershocks. Details of all focal mechanisms are given in the supplementary material (see Additional le 2). The number of polarity data ranged from 8 to 80, and its average was about 30. In the 591 focal mechanisms among 818, the P wave rst motion polarities more than 20 were used to determine the focal mechanisms. The nodal plane uncertainty ranged from 8° to 55°, and the average uncertainty for the 1636 (= 818 × 2) nodal planes of the 818 focal mechanisms was 27°. We evaluated the quality of the determined focal mechanisms as A, B, C or D based on its estimation accuracy according to Hardebeck and Sherer (2002). Quality A and D solutions have the highest and lowest levels of quality, respectively. The number of focal mechanisms of Qualities A, B, C, and D were 314, 231, 130, 143, respectively. In this study, we use the 545 focal mechanisms of Qualities A and B in the following analyses (Fig. 4). The nodal plane uncertainty of these mechanisms ranges from 8° to 42°, and the average for the 1090 (= 545 × 2) nodal planes of the 545 focal mechanisms is 22°.
The averages of the azimuths of the P-and T-axes are N78° ± 33°E and N3° ± 52°W, respectively, for all 545 focal mechanisms of Quality A and B (Fig. 5). The averages of the dips of the P-and T-axes are 25° ± 16° and 46° ± 20°, respectively. The T-axes have a larger dip angle than the P-axes. These variations of the P-and T-axes come from the variations in the focal mechanisms. Triangle diagrams (Frohlich 2001) show the distribution of focal mechanisms based on the dip angles of the P-, T-, and N-axes (Fig. 6). In the study area, most focal mechanisms are classi ed into reverse fault, strike-slip fault, and other type's earthquakes, and few focal mechanisms are classi ed into normal fault earthquakes.

Orientation of principal stresses
We obtained the stress parameters at 44 nodes in the aftershock area of the 2018 Hokkaido Eastern Iburi earthquake. The calculated values at each node are given in the supplementary material (see Additional le 3). The mean and the standard deviation of the parameters at the 44 nodes were calculated (Table 1). We found that the axis of σ 1 is oriented to ENE-WSW and the axis of σ 1 is close to be horizontal. We also found that the axis of σ 3 is close to be vertical. Therefore, the results strongly suggest that the reverse fault type stress eld is dominant, and the near-horizontal compressional stress is acting in the ENE-WSW direction in the aftershock area of the 2018 Hokkaido Eastern Iburi earthquake. According to JMA (2018a), the orientation of P-and T-axes of the CMT solution of the main shock are N67°E and N86°W, respectively, and the dip of P-and T-axes are 17° and 71°, respectively. Therefore, the stress eld obtained in this study is consistent with the CMT solution of the main shock. To investigate the stress eld in detail within the aftershock area, the deviation from the mean azimuth of σ 1 -axis is shown in Fig. 7. The mean azimuth of σ 1 -axis is N73°E and there are no remarkable spatial changes in the deviation from the mean azimuth. On the other hand, there are remarkable spatial changes in the deviation from the mean dip of σ 1 -axis. The mean dip of σ 1 -axis is 17°. The dip is deviated by approximately 10° in the direction that dip becomes horizontal in the portion shallower than 20 km and the dip is deviated by approximately 10° in the direction that dip becomes vertical in the portion deeper than 20 km. A similar tendency is seen in the dip of σ 3 -axis. In the case of the azimuth of σ 3 -axis, the spatial pattern changes at a depth of 30 km rather than 20 km.

Depth dependence of the stress ratio R
The stress ratio R is also not uniform over the aftershock area (Fig. 8). There seems to be depth dependence: R decreases systematically from the shallow to the deep portions. To con rm the depth dependency, we estimated a best-tted curve of R as a function of the depth. The nodes in the depth direction are located from a depth of 8.2 km to 45.1 km with an interval of 4.1 km. We calculated the average value of R at each depth ( Fig. 9). By tting the polynomials of 0th to 3rd to the average value of R on Table 2, AIC was calculated, and the optimal order of the polynomial was determined: where R is the stress ratio averaged at each depth and z is the depth in km. The result of the polynomial tting to R was shown on Table 2. AIC is the smallest when m = 2, therefore the depth dependence of R is not linear but quadratic.   The maximum shear stress is de ned as τ max = (σ 1 -σ 3 ) / 2. Based on the depth dependency of R, we estimated τ max as a function of depth with assumptions as follows: (1) R is given by a quadratic polynomial of depth as described above, (2)  we assumed σ 2 = 1.01 σ 3 . As a result of the calculation, τ max monotonically increases up to a depth of 27 km, reaches a maximum value, and then decreases below 27 km (Fig. 9). Note that the important point is not the absolute value of τ max , but the change pattern of increase/maximum value/decrease. The absolute value depends on how you suppose the relationship between σ 2 and σ 3 .

Discussion
Reverse-faulting stress eld in the western boundary of the HCZ In this study, the state of stress was revealed in the aftershock area of the 2018 Hokkaido Eastern Iburi earthquake (M JMA = 6.7) which is in the western boundary of the Hidaka Collision Zone (HCZ). The state of stress revealed by a stress inversion analysis of the aftershocks showed that the dominant stress eld is the reverse fault type, the σ 1 -axis is in the direction of ENE-WSW, i.e., N73°E, and the σ 1 -axis is close to be horizontal. The direction of ENE-WSW is clearly different from the convergence direction of the PA plate. Therefore, a model that the compressional stress eld due to the collision extends to the western boundary of the HCZ is strongly supported by not only the CMT solution of the 2018 main shock but also the focal mechanisms of its aftershocks. and 20°, respectively. Although the strike of σ 1 -axis obtained by Kita et al (2012) is almost same as that obtained in this study, the dip direction is different each other. This difference might suggest the stress eld is not uniform in the HCZ.

Possible causes of stress inhomogeneity
In this study we found that the stress ration R decreases as a quadratic polynomial of depth from 8 to 45 km and this change is due to the maximum shear stress τ max that increases from 8 to 27 km and decreases from 27 to 45 km. There are three possible causes of the stress inhomogeneity. The rst one is change in materials around the Moho discontinuity. The depth of the Moho was estimated to be approximately 30 km in the aftershock area of the 2018 Hokkaido Eastern Iburi earthquake (Yoshii 1972;Matsubara et al 2017). The change of τ max at a depth of 27 km from increase to decrease seems to be related to the Moho. However, mantle materials have a higher rigidity than crustal materials in general, and thus it is unlikely that the shear stress will suddenly decrease in the mantle. The second one is the brittle-ductile transition. The temperature was estimated to be approximately 400 °C at a depth of 30 km in the aftershock area and this temperature is lower than the solidus temperature (Nishida and Hashimoto 2007).
The third one is the coseismic slip during the main shock. Many authors pointed out that aftershocks concentrate not within the large coseismic slip area but in its surrounding area (e.g. Mendoza

Conclusions
We deployed temporary seismic stations immediately after the main shock of the 2018 Hokkaido Eastern Iburi earthquake (M JMA = 6.7). The dense seismic stations enabled us to determine focal mechanism solutions accurately by using the rst motion polarity of P wave. A stress inversion method was applied to the focal mechanism solutions to investigate the state of stress in the aftershock area. Looking at the stress eld in the aftershock area in detail, the orientation of σ 1 -and σ 3 -axes seemed to slightly depend on the depth. Moreover, what is interesting is the depth dependence of the stress ratio R. We presented a model that the R value is a quadratic polynomial of depth, indicating the change in shear stress that has a maximum around 27 km in depth. Although the coseismic slip seems to best explain the depth dependence of the shear stress, this model is a hypothesis to examine in future works.

Declarations
Availability of data and materials The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.