Focal mechanisms and the stress field in the aftershock area of the 2018 Hokkaido Eastern Iburi earthquake (M_JMA = 6.7)

Abstract

The tectonic stress field 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, 894 focal mechanism solutions were determined using the first-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 trend of P- and T-axes is 78° ± 33° and 352° ± 51°, respectively, and (3) the average plunge of P- and T-axes is 25° ± 16° and 44° ± 20°, respectively: the P-axis is close to be horizontal and the T-axis is more vertical than the average of the P-axes. We applied a stress inversion method to the focal mechanism solutions to estimate a stress field in the aftershock area. As a result, we found that the reverse fault type stress field is dominant in the aftershock area. An axis of the maximum principal stress (σ₁) has the trend of 72° ± 7° and the dipping eastward of 19° ± 4° and an axis of the intermediate principal stress (σ₂) has the trend of 131° ± 73° and the dipping southward of 10° ± 9°, indicating that both of σ₁- and σ₂-axes are close to be horizontal. An axis of the minimum principal stress (σ₃) has the dipping westward of 67° ± 6° that is close to be vertical. The results strongly suggest that the reverse-fault-type stress field is predominant as an average over the aftershock area which is in the western boundary of the Hidaka Collision Zone. The average of the stress ratio R = (σ₁ − σ₂)/(σ₁ − σ₃) is 0.61 ± 0.13 in the whole aftershock area. Although not statistically significant, we suggest that 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 Pacific (PA) plate is moving toward N63° W with a speed of 8.2 cm/year (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 2002, 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 field 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, 1986, 1996; Seno 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/year 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 field 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 al 1996; Moriya et al 1998; Iwasaki et al 2004; Shiina et al 2018).

Fig. 1

Map showing the Hokkaido corner. Inset shows the plate boundaries (Bird 2003) and the study area in rectangle. The region in and around the Hidaka Mountain Range (HMR) is the Hidaka Collision Zone (HCZ). Two broken lines show the eastern and the western boundaries of the HCZ. Closed triangles indicate active volcanoes. A red beach ball labeled as 2018 is the centroid moment tensor (CMT) solution of the 2018 Hokkaido Eastern Iburi earthquake (M6.7) (JMA 2018a). A red beach ball labeled as 1982 is the focal mechanism solution of the 1982 Urakawa-oki earthquake (M7.1) by using the first motion polarity of P wave (Moriya et al 1983). A red beach ball labeled as 1970 is the focal mechanism solution of the 1970 Hidaka earthquake (M6.7) by using the first motion polarity of P wave (Kita et al 2012). PA: Pacific plate, PH: Philippine Sea plate, EU: Eurasian plate, NA: North American plate, and KT: Kurile Trench

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 a complicated stress field and produce earthquakes with a variety of the focal mechanisms.

Recently, two large earthquakes occurred in the crust of the HCZ: the 1970 Hidaka earthquake (M_JMA = 6.7) (Motoya and Kitagamae 1971; Moriya 1972) and the 1982 Urakawa-oki earthquake (M_JMA = 7.1) (Moriya et al 1983). The focal mechanism solutions of the two earthquakes are similar: the reverse faulting with the P-axis in the direction of NE–SW (Kita et al 2012). Both earthquakes were in the central part of the HCZ and thus the compressional stress filed in the direction of NE–SW is dominant at least in the central part of the HCZ. This idea was supported by Terakawa and Matsu’ura (2010), founding that the reverse-fault-type stress field is dominant in and around the HCZ and the maximum principal stress (σ₁) is oriented to nearly NE–SW. Kita et al (2012) applied a stress inversion analysis to small earthquakes in the central part of the HCZ and they found that the σ₁-axis is oriented to nearly NE–SW.

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 first 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 P-axis 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 field due to the collision extends to the western boundary of the HCZ (Terakawa and Matsu’ura 2010; Hua et al 2019).

The main shock of the 2018 Hokkaido Eastern Iburi earthquake was followed by many 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 was 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 field 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; Additional file 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 telemetered systems and 22 portable offline 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 first motion polarities of P wave were read manually by a well-trained person.

Fig. 2

a Seismic stations used in this study. Temporary and permanent stations are shown in crosses and in squares, respectively. A red star indicates the epicenter of the main shock of the 2018 Hokkaido Eastern Iburi earthquake. A rectangle indicates an area shown in Figs. 3, 4, 5. Thin solid lines are active faults (Research Group for Active Faults in Japan 1991). b P wave velocity structures. A red line is based on Kasahara et al (1994) and a black line is based on Iwasaki et al (2004) and Katsumata et al (2006)

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 $V_{P} /V_{S} \, = \,\sqrt 3$, where V_P and V_S are the P and S wave velocities, respectively. We located 1785 earthquakes in the study area (42.5–42.9° N, 141.8–142.2° E), observed from 2018-09-06 03:00 to 2018-10-31 23:59 with depths shallower than 50 km and the magnitude ranging from 0.2 to 5.9 (Fig. 3). The maximum amplitude on the vertical component was measured and the magnitude was calculated by using an empirical equation (Watanabe 1971).

Fig. 3

Hypocenter distribution of aftershocks of the 2018 Hokkaido Eastern Iburi earthquake. a Epicenter distribution of aftershocks. A red star indicates the epicenter of the main shock. Thin solid lines are active faults (Research Group for Active Faults in Japan 1991). b Vertical cross section in East–West. A red star indicates the hypocenter of the main shock. c Vertical cross section in North–South. A red star indicates the hypocenter of the main shock

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 first 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) and the velocity in the crust deeper than 10 km is based on a seismic tomography analysis (Katsumata et al 2006). The second one is the same velocity structure as used in the hypocenter determination. The major difference between the two structures is the presence or absence of low-velocity sedimentary layers in the shallow part. No amplitude data were used.

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 (trend tr and plunge pl angle) of the axes of three principal stresses and the stress ratio. The trend tr is an azimuthal angle: tr = 0°, 90°, 180°, and 270° indicate the north, the east, the south, and the west, respectively. The plunge is measured from the horizontal: pl = 0° and 90° indicate horizontal and vertical axes, respectively. The principal stresses are the maximum principal stress (σ₁), the intermediate principal stress (σ₂), and the minimum principal stress (σ₃) and the stress ratio R is defined as R = (σ₁ − σ₂)/(σ₁ − σ₃), indicating the relative magnitude of the principal stresses and ranging from 0 to 1. We determined the stress field 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 sufficiently fine 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 fixed 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. The stress parameters were calculated at nodes with at least 8 focal mechanisms. 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.

Results

Focal mechanisms

We determined 894 focal mechanisms from 1785 aftershocks. Details of all focal mechanisms are given in the supplementary material (see Additional file 2). The number of polarity data ranged from 8 to 79, and its average was about 30. There were 589 focal mechanisms with more than 20 first motion polarities. The nodal plane uncertainty ranged from 5° to 70°, and the average uncertainty for the 1788 (= 894 × 2) nodal planes of the 894 focal mechanisms was 29°. 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 234, 271, 192, and 197, respectively. In this study, we use the 505 focal mechanisms of Qualities A and B in the following analyses (Fig. 4). The nodal plane uncertainty of these mechanisms ranges from 5° to 41°, and the average for the 1010 (= 505 × 2) nodal planes of the 505 focal mechanisms is 22°.

Fig. 4

Focal mechanism solutions with Quality A and B determined in this study. The depth of hypocenter is a 0–20 km, b 20–30 km, c 30–31 km, d 31–32 km, e 32–33 km, f 33–34 km, g 34–35 km, h 35–36 km, and i 36–45 km. All diagrams are equal area projections on the lower hemisphere of the focal sphere. Colored areas show compressional quadrants

The averages of the trends of the P- and T-axes are tr = 78° ± 33° and 352° ± 51°, respectively, for all 505 focal mechanisms of Quality A and B (Fig. 5). The averages of the plunges of the P- and T-axes are pl = 25° ± 16° and 44° ± 20°, respectively. The T-axes have a larger plunge 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 plunge of the P-, T-, and B-axes (Fig. 6). In the study area, most focal mechanisms are classified into reverse fault, strike-slip fault, and other type’s earthquakes, and few focal mechanisms are classified into normal fault earthquakes.

Fig. 5

Orientation of a P- and d T-axes of the focal mechanisms with Quality A and B shown in Fig. 4. Red lines in a and blue lines in d indicate the trend of the P- and T-axes, respectively. Histograms of the trend of b P- and e T-axes. Histograms of the plunge of c P- and f T-axes. The plunge of 0° and 90° indicate horizontal and vertical axes, respectively

Fig. 6

Triangle diagram displaying the distribution of the focal mechanisms with Quality A and B shown in Fig. 4. Curved lines are boundaries with T, B, and P axes of 40°, 30°, and 30° from vertical, respectively (Frohlich 2001)

Orientation of principal stresses

We performed the stress inversion using the focal mechanisms in Fig. 4 and obtained the stress parameters at 32 nodes in the aftershock area of the 2018 Hokkaido Eastern Iburi earthquake. However, some focal mechanisms in Fig. 4 were located far from all nodes and they were not used for the inversion. The calculated values at each node are given in the supplementary material (see Additional file 3). The mean and the standard deviation of the parameters at the 32 nodes were calculated (Fig. 7 and Table 1). We found that the axis of σ₁ is oriented to ENE–WSW and the axis of σ₁ is close to be horizontal or tilting down to the eastward. The fault plane of the main shock dips approximately 70° eastward (Kobayashi et al. 2019a,b; Guo et al. 2019). The pl = 19° of σ₁ may promote the reverse faulting slip of the main shock. We also found that the axis of σ₃ is close to be vertical. Therefore, the results strongly suggest that the reverse fault type stress field 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 trend of P- and T-axes of the CMT solution of the main shock are tr = 67° and 274°, respectively, and the plunge of P- and T-axes are pl = 17° and 71°, respectively. Therefore, the stress field obtained in this study is consistent with the CMT solution of the main shock. The vertical cross sections are shown in Fig. 8. No stress parameter was obtained in the aftershock area shallower than 20 km except for one node. There is no remarkable systematic spatial variation in trend and plunge of the σ₁- and σ₃-axes in the aftershock area deeper than 20 km. This observation suggests that at least the orientation of the principal stresses is uniform in the aftershock area.

Fig. 7

Stress parameters at nodes shown in Fig. 8. a The orientation of principal stresses: red circles, crosses, and blue circles indicate σ₁, σ₂, and σ₃ at each node, respectively. b Histogram of the stress ratio R = (σ₁—σ₂)/(σ₁—σ₃)

Table 1 Principal stresses averaged over the aftershock area

Fig. 8

Spatial pattern of the trend and the plunge. a The trend of σ₁-axis, b the plunge of σ₁-axis, c the trend of σ₃-axis, and d the plunge of σ₃-axis. The coseismic slip was indicated by contours with 0.6 and 1.2 m (Asano and Iwata 2019)

Depth dependence of the stress ratio R

When the 32 nodes were averaged, the stress ratio R was 0.61 ± 0.13 (Table 1). Ohtani and Imanishi (2019) conducted a stress inversion using 27 focal mechanisms in the aftershock area and obtained the stress ratio φ = 1.0–R = 0.57. Since they did not describe the exact value of φ in the text, the approximate values were read from a histogram in Fig. 2b of Ohtani and Imanishi (2019). The value φ = 0.57, that is, R = 0.43 is almost in the range of 1σ (= 0.13) and thus we cannot say statistically that the stress ratio R = 0.43 obtained by Ohtani and Imanishi (2019) is different from R = 0.61 obtained in this study. Consequently, from the point of view of the average of the entire aftershock area, there is no significant difference between the result of Ohtani and Imanishi (2019) and ours.

The estimated value of R at each node has a very large error (see Additional file 3). The difference in the values between the nodes falls within the error range, and it cannot be said that there is a statistically significant difference. Therefore, the following judgment is reasonable: it is futile to discuss further the spatial pattern of stress ratios. Although we know that there is no statistical significance, we dare to propose a hypothesis that R may have depth dependence in this study, as seen in the averaged values of R.

There seems to be depth dependence: R decreases systematically from the shallow to the deep portions (Fig. 9). To express the depth dependency quantitatively, we estimated a best-fitted curve of R as a function of the depth. The nodes in the depth direction are located from a depth of 16.4–45.1 km with an interval of 4.1 km. We calculated the average value of R at each depth (Fig. 10a). By fitting the polynomials of 0th to 3rd to the average value of R, AIC was calculated, and the optimal order of the polynomial was determined:

$$R\, = \,\mathop \sum \limits_{n\, = \,0}^{m} a_{n} z^{n} \;\left( {m = \,0, 1, 2, 3} \right)$$

(1)

Fig. 9

Distribution of the stress ratio R in vertical cross section. The coseismic slip was indicated by contours with 0.6 and 1.2 m (Asano and Iwata 2019)

Fig. 10

Depth dependency of the stress ratio R. a Black dots are the R-value at each node and red crosses are the averaged values at each depth. Horizontal thin lines indicate the estimation error. A curved line is the best fitted line to the red crosses. b The maximum shear stress τ_max = (σ₁−σ₃)/2 as a function of depth. c Histogram of the depth of aftershocks shown in Fig. 3

where R is the stress ratio averaged at each depth and z is the depth in km. The result of the polynomial fitting 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.

Table 2 Polynomial fitting to the stress ratio R

The maximum shear stress is defined as τ_max = (σ₁ − σ₃)/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) the minimum principal stress σ₃ is equal to the lithostatic overburden pressure minus hydrostatic pressure, σ₃(z) = 16.7 z (MPa) at z km depth (e.g. Aochi and Ulrich 2015; Ando and Kaneko 2018; Hisakawa et al 2020), and (3) σ₂ ≈ σ₃ (Hisakawa et al 2020) in this case we assumed σ₂ = 1.01 σ₃. As a result of a simple arithmetic calculation, τ_max monotonically increases up to a depth of 32 km, reaches a maximum value, and then decreases below 32 km (Fig. 10b). Interestingly, the aftershock activity concentrates around the depth of 32 km (Fig. 10c). 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 σ₂ and σ₃. For example, if you assumed σ₂ = 1.1 σ₃, τ_max is 10 times larger.

Discussion

Reverse-faulting stress field 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 field is the reverse fault type, the σ₁-axis is in the direction of ENE-WSW, i.e., tr = 72°, and the σ₁-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 field 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.

Terakawa and Matsu’ura (2010) insisted that the stress field of reverse faulting type is dominant in the aftershock area of the 2018 Hokkaido Eastern Iburi earthquake: the axis of σ₁ is close to be horizontal with almost EW direction and the axis of σ₃ is near-vertical. The results obtained in this study are consistent with those obtained by Terakawa and Matsu’ura (2010). Kita et al (2012). applied a stress inversion analysis to focal mechanisms of small earthquakes occurred in and around an area 100 km east of the 2018 main shock. As a result, they found that the trend and the plunge of σ₁-axis are 224° and 20°, respectively. The trend of σ₁-axis obtained by Kita et al (2012) is different from that obtained in this study. This difference might suggest the stress field is not uniform in the HCZ.

A possible cause of the aftershock concentration in the deeper zone than 30 km

Many authors pointed out that aftershocks concentrate not within the large coseismic slip area but in its surrounding area (e.g. Mendoza and Hartzell 1988; Beroza and Zoback 1993; Das and Henry 2003; Hsu et al. 2006; Woessner et al. 2006; Perfettini et al. 2010; Asano et al. 2011; Kato and Igarashi 2012). These observations strongly suggest that aftershocks are induced by the local accumulation of the shear stress due to the rupture of main shock.

The aftershock activity of the 2018 Hokkaido Eastern Iburi earthquake concentrates between 30 and 40 km in depth. The depth of D90 is approximately 25 km in this area (Omuralieva et al 2012). The depth of D90 is defined as the depth above which 90% of the earthquakes occur. The aftershock concentration zone is much deeper than the depth of D90. In this study we suggested that the stress ratio R decreases as a quadratic polynomial of depth from 16 to 45 km and this change is due to the maximum shear stress τ_max that increases from 16 to 32 km and decreases from 32 to 45 km. The aftershock activity has concentrated in the zone that the τ_max becomes the maximum. Some authors have reported the coseismic slip distribution and showed that the seismic fault ruptured by the main shock did not expand deeper than 30 km (Kobayashi et al. 2019a,b; Asano and Iwata 2019). Asano and Iwata (2019) analyzed the strong-motion data by using a kinematic waveform inversion method and found that the peak slip of 1.7 m was located at a depth of about 26 km, southwest of the epicenter.

A key to explain these observations is the brittle–ductile transition zone. The depth to the Moho is approximately 30 km in the aftershock area (Yoshii 1972; Matsubara et al 2017) and the temperature at the Moho is estimated to be 400—600 °C (Fujiwara 1984; Nishida and Hashimoto 2007). The deeper part of the aftershock area than 30 km is likely to be the brittle–ductile transition zone. Although a long-term stress of thousands of years does not accumulate in the brittle-ductile transition zone, a short-term stress of several years accumulates. Therefore, we suggest that the rupture of the main shock was limited shallower than 30 km, due to this rupture, a transient shear stress was loaded in the mantle deeper than 30 km, and the aftershock activity was induced there. Hisakawa et al (2020) conducted a dynamic rupture simulation of the2018 Hokkaido Eastern Iburi earthquake and showed that the shear stress rise of 10–15 MPa is possible theoretically near the boundary on the deep side of the large cosismic slip area.

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 first 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. At deeper than 20 km, the orientation of σ₁- and σ₃-axes seemed to be uniform. Moreover, what is interesting is the depth dependence of the stress ratio R. Although the statistical significance is very poor, we dare to present a model that the R value is a quadratic polynomial of depth, indicating the change in shear stress that has a maximum around 32 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.