Dynamic rupture simulation of 2018, Hokkaido Eastern Iburi earthquake: role of non-planar geometry

The 2018, Hokkaido Eastern Iburi, Japan, earthquake is an event characterized by complexity of the rupture process and slip pattern, which may involve both reverse and strike-slip motion depending on the locations on the fault surface. We perform dynamic rupture simulations based on simple physical laws, conditions for stressing and fault friction, and the non-planar fault geometry constrained by the aftershock observation. The complex fault geometry is numerically treated by the boundary integral equation method accelerated by the fast domain portioning method. The fault geometry is characterized primarily by the combination of six fault planes. As a result, we are able to explain several observed features of the event, including the spatial variation of the final fault slip and rupture velocity, which are inferred from the kinematic slip inversion. We also succeed in refining the constraint of the regional stress field in the focal area based on the simulation. Our results show that the overall patterns of the complex rupture event can be reproduced by a relatively simple model of the regional stress and the fault friction, if the geometrical complexity of the fault is properly taken into account.


Introduction
The 2018, Hokkaido Eastern Iburi earthquake occurred on September 6 (JST) with an estimated local magnitude ( M JMA ) of 6.7. The event occurred in so-called the Hokkaido corner region, where the northern Japan island arc is considered to transit to the Kurile island arc (Kita et al. 2012;Kimura 1994). The Kurile island arc is moving in the southwest direction and is colliding with the Japan island arc (Kita et al. 2012;Kimura 1994), which causes the overall tendency of the stress field in the region to be northeast-southwest compression (Fig. 1 inset).
The centroid moment tensor (CMT) solution determined by Japan Meteorological Agency (JMA) (https :// www.data.jma.go.jp/svd/eqev/data/mech/cmt/fig/cmt20 18090 60307 59.html) showed that the event ruptured east dipping reverse faults with some non-double-couple component (Fig. 1a). The first-motion solution determined by JMA contained more strike-slip components than the CMT solution. This indicates that the rake angle has changed during the rupture process. As shown in Fig. 1a, the bent aftershock distribution (National Research Institute for Earth Science and Disaster Resilience (NIED) 2019a; Matsubara et al. 2019) clearly indicates that the fault is not a single plane, which is a clue for understanding the origin of complexity in the rupture process.
In this study, we first aim to explain the physical reasons for this complex rupture process. While previous studies provided insights into the kinematic feature of this earthquake (Asano and Iwata 2019), physical conditions such as stress field underlying the event are poorly understood. Therefore, we conducted dynamic rupture simulations by giving minimal assumptions constrained with observed facts (i.e., the aftershock distribution, the stress tensor inversion from CMT solutions, the principal stress σ v proportional to depth z (Brudy et al. 1997), the first motion polarity) as much as we could, rather than introducing many tuning parameters, which makes the model validation difficult. In this sense, we do not use the slip inversion results to determine the key parameters of the stress drop and frictional strengths and the fault geometry.
The second purpose of this study is to re-constrain the local stress field. The stress field has been inverted from the focal mechanisms of the middle-sized earthquakes in the region (Terakawa and Matsu'ura 2010), but the result is insufficient to resolve the specific extent of the focal area. Therefore, we test many sets of stress parameters to determine which ones produce a model that best fits the observations. We also discuss the robustness of those stresses and their fit. We make use of the fully dynamic 3-D rupture simulation based on the spatiotemporal boundary integral equation method (ST-BIEM) (e.g., Aochi and Fukuyama 2002;Tada 2006;Ando and Kaneko 2018). While the 3-D non-planar fault dynamics can be analyzed by other numerical methods such as the finite element method (e.g., Oglesby and Mai, 2012), the spectral element method (e.g., Galvez et al. 2014) and the discontinuous Galerkin method (e.g., Pelties et al. 2012), the strong point of this boundary-based method is its high numerical accuracy and the easiness of handling non-planar complex fault geometry including branching by meshing without considering the surrounding medium. The recent development  (Kita et al. 2012), respectively. a The map view of aftershock distribution. The red star indicates the observed mainshock epicenter. The red lines show the cross-sectional projection of the assumed fault at a depth of 34 km. The orange line shows the branch fault (Segment N) added in Case B. b The cross-sectional views for the sections i-v; the red and orange lines show the cross-sectional projection of the assumed fault surfaces on the middle of each section, respectively. The red star indicates the locations of the observed hypocenter and the assumed hypocenter in Case B, and the yellow star indicates the assumed nucleation location in Case A without the branch fault of the fast domain partitioning BIEM (FDP-BIEM) by Ando (2016) has allowed us to conduct the analysis at the reduced cost of O N 2 for the given number of elements N.

Numerical method
The dynamic rupture process can be solved as a forward problem by assuming the fault geometry, the stress field, and the physical laws such as the friction law and the elastic response of the medium. In FDP-BIEM, we solve the discretized boundary integral equations (BIEs) obtained from representation theorem in real space and time domain, representing the relationship between the slip history and the current traction in the form of where T i,n p denotes the pth component of the traction on the ith fault element at nth time step, T i,init p the initial traction, D j,m q the qth component of the slip rate on the jth element of fault at the past mth time step, K i,j,n−m pq the integration kernel representing the elastodynamic responses due to the past slip on each fault element, and the constants µ and β are rigidity and S-wave speed, respectively. We do not restrict the slip direction (q = 1 for the strike-slip or 2 for the dip-slip components) and the slip direction can be temporal variable determined by the transient direction of the traction. The triangular boundary elements (Tada 2006) are employed. The fast domain partitioning method is used for the efficient summation and the compact memory storage of the kernel on the right-hand side of (1).
The boundary condition to constrain D or T for each element is given by where µ i,n fric denote the frictional coefficient on the ith fault element at nth time step. We obtain the value of µ fric to follow the linear slip-weakening friction law with characteristic slip distance D c where µ s and µ d denote static frictional coefficient and dynamic frictional coefficient, and S i,n is the slip vector calculated by S i,n = S i,n−1 + D i,n t , respectively. t (1) is the time step size. By considering (1)-(5), the equation system becomes closed and can be solved numerically. We assume no opening mode of fault displacement, which gives D i,n 3 = 0. We employed fault models discretized by triangles with an average side length of 200 m. We used Gmsh (http:// www.geuz.org/gmsh) for meshing. The fault models consist of about 35,000 triangular fault elements, and we calculated over 3000 time steps for each rupture event. The calculation for this case required about 4 TB of memory. The total calculation time in this case was about 20 min with the 512 CPUs (34,816 cores).

Model and observational constraints
It is known that the fault geometry (e.g., King 1986;Yoshida et al. 1996) together with the regional stress field (e.g., Aochi and Fukuyama 2002;Ando et al. 2017;Ando and Kaneko 2018) affect the rupture process to a large extent and we take into account the effects of the nonplanar fault surfaces and the regional stress field in our model. Therefore, to validate our model and make it testable with comparing to the slip inversion result (Asano and Iwata 2019), we constrain these parameters based on observations independent from slip inversion results.

Fault geometry model
We assume a homogeneous isotropic elastic half-space ( µ = 40 GPa, β = 4 km/s (Matsubara et al. 2017)), and place a fault. We constrained 3-D fault geometry based on the hypocentral distribution of the aftershocks within 3 days after the mainshock (NIED 2019a; Matsubara et al. 2019), relocated with the hypoDD method (Waldhauser and Ellsworth 2000) derived from High-sensitivity seismograph network (Hi-net) operated by NIED (Okada et al. 2004;Obara et al. 2005;NIED 2019b). The fault geometry is visually determined to describe the aftershock distribution for the five sections from I to V as shown in Fig. 1b.
In the section i, the alignment of the aftershocks is identified in the direction of N 12° E on the map view (Fig. 1a). In the cross-sectional view (Fig. 1b, section i), the aftershocks appear to be localized into a plane with the dip angle of 64° above the depth of about 34 km and a plane of 68° below it. We set no fault in the deeper part of this section (depth of > 43 km) because the distribution of the aftershock in the area looked scattered, not localized into a plane.
In the section ii, the southern part of the area has almost the same tendency as section i, but the aftershocks are aligned in the direction of strike N 13° W, as represented by the relatively steeper and shallower red lines located above about 34 km depth, respectively. The much steeper plane at the depth is also identified. In section iii, the fault strike is almost the same as the northern part of section ii. The assumed fault dip angles are 58° and 75°, respectively, for the shallower and deeper than the bend at 34 km depth (Fig. 1b, section iv). Between the sections iii and iv, we identify the horizontal bend of the aftershock alignment, where the strike of the alignment in section iv and v is almost the same as that in the section i.
Note that the fault surface is extended to the shallower part (until about 20 km depth) where the aftershocks are less crowded because it can be assumed that there occurred less aftershocks in the part of large slip, and released large stress even though the fault exists (Discussed later on this paper). Overall, we assume the continuous fault surface consisting of the six planes with the bends at the depth of 34 km and the longitudes of 42° 39′ N and 42° 43′ N as shown in Fig. 2. The top of the fault is at 20 km depth and the bottom of it is at 43 km depth.
In addition to the abovementioned main fault, we also consider a branch fault called Segment N ( Fig. 1a, b, section iii), which presumably hosted the actual nucleation of the mainshock. As shown in Fig. 3, Segment N is assumed to be on the western side of the main fault connected to Segment 3 and Segment 4, the size of which is about 5 km × 5 km and comparable to what is suggested by an earlier report discussing the first motion polarity (Katsumata et al. 2019). The angles of the strike (N 74° W) and dip (48°) of Segment N are determined based on the first motion reported by JMA (Fig. 1a). On Segment N, the nucleation is assumed at the junction connected to Segment 4 at a depth of 36 km as corresponding to the observed hypocenter ( Fig. 1b, section iii). Hereafter, we call the fault model without the branch fault as Case A and that with it as Case B. The effect of the different stress parameters on the slip pattern is investigated mainly using Case A, while Case B is used to examine the detailed process of the nucleation.

Regional stress field and initial conditions
We aim to constrain the regional stress field based on observational information that can be obtained before the event and on established empirical laws as below. First, we constrained the principal stress axes v i ( i = 1, 2, 3 for the maximum, intermediate, and minimum principal stress axes, respectively). We assumed that one of them to be vertical following the Andersonian condition (Anderson 1951). To constrain the remaining two horizontal principal axes, we referred to the stress inversion result around Japan obtained by Terakawa and Matsu'ura (2010). Due to their limited spatial resolution, we roughly determine the maximum principal stress (called S1 hereafter) axis to be directed to the northeast (about N 50° E) from the average of the surrounding areas. We will examine the possible ranges of S1 direction later based on the dynamic rupture simulation.
Next, we constrained the value of the maximum, intermediate, and minimum principal stress, σ 1 , σ 2 , and σ 3 . Here, we define two ratios to show the relationship of the stress values and The result of stress inversion implies that the region is located close to the transition zone between the stress field of strike-slip events and that of reverse fault events, indicating that σ 2 ≈ σ 3 and giving ϕ ≈ 0 in our definition. We assume that the minimum principal stress (S3) is in the vertical direction and its magnitude σ 3 is determined by lithostatic overburden pressure minus hydrostatic pressure, σ 3 (z) = 16.7z (MPa) at z km depth as assumed in the previous study (e.g., Aochi and Ulrich 2015;Ando and Kaneko 2018). This assumption takes account of the increasing absolute stresses observed by Brudy et al. (1997) and accounts for increasing stress drop with depth. The ratio ϕ ′ cannot be determined from the stress inversion, so the parameter is set for the potential value of the stress drop τ 0 t − µ d τ 0 n (called potential stress drop hereafter) to be about 10 MPa at the location of the hypocenter on the main fault, where τ 0 n and τ 0 t denote the initial values of the normal traction and the maximum shear traction, respectively. The value of 10 MPa is chosen as a representative value, not contradicting the seismologically established values (e.g., Kanamori and Anderson 1975), which also reproduces the amount of slip in this event as shown later on this paper. By fixing the value of the stress drop at the location of the nucleation (hypocenter), the nucleation size becomes the same among the cases of the different stress parameters. We can then simply extract the effects of the interplay between the regional stress and the fault geometry from the comparison of the obtained slip distributions and rupture patterns. Note that while we assume the homogeneous regional stress field, the traction on the fault model is heterogeneous reflecting the spatial variation of the fault geometry. The projection of the stress σ to the traction on the fault surface is conducted by calculating inner products where e i p denotes an unit vector directed to pth direction at ith element, and the variation of e i p causes heterogeneous traction or potential stress drop. Although there may be some heterogeneity in the regional stress field, it was difficult to find the higher resolution and robust results of stress inversion conducted before the current event in a literature as far as we searched. Therefore, the homogeneous stress field is assumed as the robust assumption, and we focus on the effect of heterogeneous traction due to the fault geometry.
As the initial condition, we place a circular nucleation patch imposing the static slip compatible to the constant stress drop there with the critical radius of which is determined by the balance of the fracture surface energy and the energy release rate (Andrews 1976) where τ n and τ t denote the normal traction and the maximum shear traction, respectively. This nucleation procedure allows the rupture to be propagated spontaneously. In Case A without Segment N, we assumed the nucleation at the bending point between Segment 3 and Segment 4, where the potential stress drop on the main fault is the largest on Segment 3 (Fig. 4). In Case B, the nucleation is assumed on Segment N to be deeper than the case of Case A (see Fig. 3).
To examine the relationship between the stress field and the characteristics of the rupture process, we test 9 sets of the stress parameters (called Models 1-9 as in Fig. 4) considering the ambiguities of the observational constraints. For S1 axis, we consider N 45° E, N 48° E, and N 51° E, and the stress ratio ϕ = 0.01, 0.1, and 0.2 . Since we focus on the stress field and the rupture process, we set homogeneous frictional parameters µ s = 0.25 , µ d = 0.2 , and D c = 0.1 m, and we decided not to change them. The value of µ s appears to be small but one may consider this value as that includes the effect of the overpressure of the pore fluid (e.g., Andrews 2002;Oglesby and Mai 2012). Note that the fracture energy is about 5 × 10 5 Jm −2 (needed to restrict the nucleation size to approximately 1/10 of the entire fault length), which is determined by reasonable values of D c , and the strength drop (µ s − µ d )τ 0 n , and the stress drop of about 10 MPa.
The strength excess µ s τ 0 t − τ 0 n and the potential stress drop for the nine cases of the stress parameters are shown in Fig. 4. As for the potential stress drop, we observe that the value on the central shallower part of the fault (Segment 3) is almost constant for all the stress models because the hypocenter is on the Segment 3 and ϕ ′ is determined for the stress drop of this segment at the depth to be about 10 MPa. The potential stress drop on the central deeper part (Segment 4) becomes negative because this segment is steep and oblique for any of the assumed stress fields. The potential stress drop on the other segments can be either positive or negative, depending on the stress model. Overall, the potential stress drop becomes smaller for more northward-directed S1 axes and larger ϕ . The strength excess shows a similar trend of increasing and decreasing depending on the stress parameters.

Rupture process
In this section, we focus on the case of ϕ = 0.1 and S1 axis N 48° E (Model 5) in Case A to illustrate the characteristic pattern of the rupture process on the main fault because it reproduces the feature of the inversion result in that it shows the largest slip on the middle of shallower part (segment 3) and small slip on the deeper part (segments 2 and 6). Figure 5 shows snapshots of the slip rate for Model 5 at every 2 s from t = 0 s to t = 12 s. The earlier process, including the rupture transfer from Segment N to the main fault, is discussed later with the analysis result of Case B.
First, the rupture propagates to shallower part (Segment 3), without breaking the deeper part, which corresponds to the Segment 4 presenting the negative potential stress drop. This process continues from t = 0 s to 2 s. Next, the rupture begins to propagate to the northern shallower part (Segment 5) and the southern shallower part (Segment 1). In this process, the rupture velocity is different on each segment. On Segment 1, the rupture propagates mainly to the south and the velocity is approximately 1.5 km/s. On Segment 3, it propagates upward and the velocity is approximately 4.0 km/s, which is close to the S-wave speed. On Segment 5, it propagates mainly to the north and the velocity is approximately 2.0 km/s. By comparing this result to the potential stress drop distribution (Fig. 4), it can be said that the rupture velocity has a positive correlation with the potential stress drop on these segments. This process continues from t = 2 s to 6 s. The rupture on Segment 3 reaches the top of the fault at t = 6 s, and the rupture on Segment 1 and Segment 5 reaches the southern end of the fault at t = 7 s and the northern end at t = 6.5 s, respectively.
After t = 6 s, the rupture begins to propagate toward the deeper part of the fault (Segments 2 and 6) at the rupture velocity below 1 km/s. It stops half-way of the fault both on Segments 2 and 6 (see Fig. 5, t = 10 s). It is noted that on the shallower portion of the fault, the rupture to the south is arrested on the southern deepest part of Segment 1 (t = 10 s). The rupture stops almost entirely on the fault at t = 12 s.
The reason why the rupture stopped half-way on Segments 2 and 6 even though the area has positive stress, drop can be explained by taking the energy balance into account. The crack with the length of L has available energy to extend the crack by dL (Andrews 1976) Note that this value is proportional to L. The absorbed energy at the end of the fault when the crack extends by dL (Andrews 1976) where τ 0 denotes the initial shear traction, τ u the upper yield shear traction ( τ u = µ s τ n ), and τ f the dynamicfrictional level of shear traction ( τ f = µ d τ n ). The rupture propagates when E a > E r . This indicates that the rupture is more likely to propagate when L is larger for the same stress drop and the same frictional strength. The potential stress drop on Segment 1, Segment 2, Segment 5, and Segment 6 are almost the same (about 3-5 MPa, Fig. 4a). Here, we suppose L eff as the minimum characteristic length of the crack or slipping area. On Segments 1 and 5, the rupture propagates mainly in the north-south direction, and L eff is almost equal to the along-strike length of the fault (about 30 km). On Segments 2 and 6, the rupture propagates downward, and L eff is reduced to the along-dip lengths of the segments (about 20 km), limited by the Segment 4 affected as the barrier. This smaller L eff on Segments 2 and 6 may have prevented the faults from breaking entirely.

Comparison with kinematic slip inversion
The rupture process has been inverted kinematically making use of the observational record of the strong ground motion (Asano and Iwata 2019). We compare the rupture process of the current models with this slip inversion result to validate the characteristic of our model by focusing on Model 5.
We compare the slip distribution at each time during the rupture. Figure 6 shows the temporal progress of the slip obtained from the inversion, and Fig. 7 shows the progress obtained from dynamic rupture simulation, respectively, as the map views. In the first phase, both results show that the slip (until 5 s in the inversion, 2 s in simulation) is limited in the vicinity of the hypocenter. Then, the slip mainly propagated bilaterally toward north and south in the shallow portion (western side in the map view) in both of them (6-10 s in the inversion, 2-8 s in the simulation). The slip inversion result shows that the rupture velocity is faster to the north than to the south as indicated in Fig. 6. Finally, the rupture propagates to the southern deeper part (10-13 s in the inversion, 8-12 s in the simulation) and is arrested.
The results have several features in common in rupture process and slip profiles in the final distribution of the slip (15 s in the inversion and 12 s in the simulation). First, the maximum slip amount is 1.7 m in the inversion and 1.9 m in the simulation, respectively. The maximum slip also takes place almost at the same part of the fault, up-dip of the hypocenter to the south, where the fault dip is shallower becoming optimal to the regional stress field. Second, there is a tendency of the slip to be larger in the shallower part than in the deeper part. Regarding the rupture propagation patterns, the rupture velocities are faster for the propagation to the north than the south commonly in the inversion and the simulation.
The slip inversion result also shows the rake angles at 121° on the shallower potion of the fault corresponding to the location where the largest slip took place (circled at   Yagi and Okuwaki (2018) also shows the similar spatial variation of the rake angles. Such spatial distribution is clearly seen in the simulation (Fig. 7, inset) as the effect of different dip angles of the shallower and deeper segments.
We can still see a difference in the deepest part (38-42 km depth) of the fault slip particularly near the southern end. While the deepest part does not rupture (arrested at a half-way) in the simulation, the inversion result appears to propagate deeper. The cause of this difference is unclear, and it could be caused by some finer scaled fault heterogeneity not considered in our model or Blue arrows indicate the outline of the rupture propagation. The red lines drawn from t = 6 to 9 s, respectively, connect the rupture fronts propagating to north and south; the larger inclination of the north indicates faster rupture velocity this is regarded within the uncertainty of the inversion as discussed in Asano and Iwata (2019).
Next, we compare a simulated moment rate function (MRF) to that derived from the inversion result. The MRF derived from the inversion and the simulation of two models (Models 5 and 7) are shown in Fig. 8. Each plotted MRF has been normalized and shifted in time to compare the shape of their peaks. The shapes of three MRFs seem to have a similar rise from 0 to 7 s. After that the values of three MRFs begin to decrease, but the decreasing rate is different; Model 7 shows the steepest decrease, the inversion the gentlest, Model 5 the intermediate. The shorter durations of the simulation than the inversion may primarily reflect the slightly shorter fault length assumed in our model. The difference of the declining rate also reflects the difference of the slip on the deeper portion of the fault (see Figs. 6 and 7), and as shown in the simulation for Model 5 (Fig. 7), the rupture on the deeper part occurring after 8 s prolongs the declining phase (Model 7 presents smaller deep slip as shown later).
The rapid increase of the MRF in the simulation can be influenced by the assumed nucleation. For the simulation, an initial crack with a radius of about 2 km (~ M w 5) is assumed, and the process for the crack to grow to that size is ignored. The location of the nucleation can be also affective. We set the nucleation up-dip from the bend in the fault where the stress drop in the model is the largest over the fault area, though the observationally determined hypocenter can be down-dip from the bend (Fig. 1b, section iii). Propagating over the fault bend or the rupture transfer from a branch fault can take longer time resulted in prolonged rise. In the real fault, there may be finer scaled heterogeneity not included in this model (Case A), which might have allowed the rupture to nucleate on the deeper part of the fault. In the next, we examine the nucleation process with Case B taking account of the branch fault that could host the nucleation.

Effect of the nucleation location: addition of a branch fault (Segment N)
To discuss the more detailed nucleation process, we examined the behavior of Case B containing the branch fault (Segment N). We conducted a dynamic rupture simulation under the stress field of Model 7, which is the optimal stress condition as we describe later on this paper. The potential stress drop in the nucleation zone for Model 7 is about 26 MPa, meaning that the nucleation is easier to occur on Segment N than Segment 3. The potential stress drop for the other stress models is shown in Table 1 and shows similar values, which are larger on Segment N than on Segment 3. This result indicates that a qualitatively similar nucleation process will be observed in the other stress models, and thus we chose Model 7 as a representative model.
As a result of the simulation, the first motion was reproduced to a large extent (Fig. 9a), although we did not specify the rake angle while the strike and dip angles are only assumed based on the JMA focal mechanism. The duration was about 2 s from the nucleation to the propagation to its junction with Segment 3, which is the nucleation location assumed in Case A. The oblique reverse fault motion of Segment N increased the shear traction and decreased the normal traction on a specific area on Segment 3. These prolonged initial processes may explain the time difference of the MRFs and the rupture process.
The distribution of the final slip on the main fault in Case B is shown to be similar to that of Case A (see Figs. 9b and 11) in that it shows the largest slip of about 2.0 m on Segment 3, and Segments 5 and 1 are respectively ruptured entirely and partly. Also, there is little slip on deeper segments. This result clearly shows that the rupture on Segment N does not have a significant effect on the final slip pattern and rupture process on the main fault. Thus, we can safely use the results of Case A to investigate the dominant rupture processes, which occurred on the main fault.   et al. Earth, Planets and Space (2020) 72:36 Implication to spatial variation of aftershock activity Aftershocks of this event are more active in the deeper part of the fault than the shallower part even though there was little large slip in the deeper part (see Figs. 1b,6,and 7). This might be interpreted as the result of the rupture arrest on the half-way of the fault and the remained stress concentration at the neutral criticality condition. Figure 10 shows the potential stress drop on the fault before the event and after the event. In the shallower part, where the slip was large, the potential stress drop decreases to almost zero or shows some negative value because a large amount of stress has been released during the event. On the other hand, on the deeper part of the fault remaining without being ruptured, the potential stress drop increased by 10-15 MPa near the rupture front (Fig. 10b). This increase could have led to occurrence of intensive aftershocks there.
Dependence and sensitivity of the rupture processes on the stress models Figure 11 shows the dependence of the final slip distribution on the different stress models, Fig. 12 the rake distribution, and Fig. 13 the CMT solutions, for the fault geometry of Case A. Table 2 summarizes the obtained characteristics of each model regarding the slip amount and rake angles. We first observe that there is little rupture on Segment 4, reflecting the negative stress drop on the fault in any model. In each model, we find that the slip contains a more strike-slip component on Segment 2 or Segment 6, than on Segment 1 or Segment 5, reflecting the steeper dip angle, and the slip on Segments 1 and 5 contains more strike-slip component than that on Segment 3 by the effect of the fault strike. If we compare the cases resulted in the similar slip distribution (c.f., Model 1 and Model 2; Model 8 and Model 9; Model 6 and Model 9), the slip contains a more strike-slip component for smaller values of ϕ and S1 axis directed more to the north. This increased strike-slip component is understandable since smaller values of ϕ decrease the minimum horizontal compressive stress (S2). The more northward the compressive direction, the lower the angle between S1 and the fault strike, which is also closer to the optimal angle for strike-slip movement.
On the other hand, the spatial distribution of the rake angle is not as simple when the slip distribution is different. When we compare the slip on Segment 3 of Models 4, 5 and 6, we observe that Model 6 contains the least strike-slip component, which seems contradictory for the result of the above similar slip distribution. This is probably because the slip on Segment 3 of Models 4 and 5 is affected by the slip on Segments 1 or 5, which contains more strike-slip component.
Next, we examine which model explains the real event best. The slip distribution derived from the inversion shows the largest slip in the shallower part of the fault. It also shows that the shallower part of the fault has been ruptured almost entirely. In the result of each model, this feature of slip distribution is well reproduced by Models 5 and 7. The moment magnitudes in these models are 6.7, which are also consistent with the kinematic results (Asano and Iwata 2019). Among Models 5 and 7, the slip distribution is better explained by Model 5 because it shows relatively larger slip on the deeper part of the fault,   as indicated from the slip inversion. The MRF is also reproduced better in Model 5 (Fig. 8), associated with better reproduced slip on the deeper part.
The observationally constrained CMT solution by F-net is also shown in Fig. 13. We can see the observationally obtained CMT is slightly better reproduced by Model 7 than Model 5 with larger reverse faulting component although the observed CMT contains a little more reverse component.
To summarize, the above parameter study constrains the possible stress state to Models 5 and 7, where the MRF and the rupture pattern on the deeper part have been reproduced better by Model 5, and the CMT solution and the rupture on the shallower part are better with Model 7. This observed tendency may indicate that the stress state is more like Model 7 in the shallower part and is more like Model 5 in the deeper part.

Conclusions
We investigate the origin of the complex pattern observed in the CMT solutions, the slip profile and the rupture process of the 2018, Hokkaido Eastern Iburi earthquake and what kind of stress field underlying the observed characteristics. We conducted a set of dynamic rupture simulation by assuming the fault geometry observationally constrained by the aftershock distribution, the stress field, and homogeneously distributed frictional coefficients. The simulation explains significant features of the mainshock, including the overall distribution of the slip amount and rake angles on the fault, nondouble-couple component of CMT, and the difference of rupture velocity depending on the rupture direction. The regional stress field in this area is also constrained more in detailed than the previous inference with the parameter study regarding the regional stress model. We show that the observed complexity of the event can be explained primarily by the effect of non-planar fault geometry with multiple bends.

Table 2 Rake angle (on Segment 3), max slip, and the ruptured segments in each model
For the max slip, the segment on which the slip occurs is also displayed. For the ruptured segment, segments put in parentheses are ruptured partly