Estimation of the source process and forward simulation of long-period ground motion of the 2018 Hokkaido Eastern Iburi, Japan, earthquake

In this study, we investigate the source rupture process of the 2018 Hokkaido Eastern Iburi earthquake in Japan (MJMA 6.7) and how the ground motion can be reproduced using available source and velocity models. First, we conduct a multiple-time-window kinematic waveform inversion using strong motion waveforms, which indicates that a large slip area located at a depth of 25–30 km in the up-dip direction from the hypocenter was caused by a rupture propagating upward 6–12 s after its initiation. Moreover, the high-seismicity area of aftershocks did not overlap with the large-slip area. Subsequently, using the obtained source model and a three-dimensional velocity structure model, we conduct a forward long-period (< 0.5 Hz) ground-motion simulation. The simulation was able to reproduce the overall ground-motion characteristics in the sedimentary layers of the Ishikari Lowland.


Background
The 2018 Hokkaido Eastern Iburi earthquake (called "the 2018 Iburi earthquake" hereafter) occurred in the eastern Iburi region of Hokkaido, Japan, at 3:08 JST on September 6, 2018 (18:08 UTC on September 5, 2018) (Fig. 1a). The Japan Meteorological Agency (JMA) magnitude (M JMA ) is 6.7. The moment magnitude estimated from moment tensor analysis of F-net of the National Research Institute for Earth Science and Disaster Resilience (NIED) (Fukuyama et al. 1998) and that obtained by the Global centroid moment tensor (GCMT) project are 6.6 and 6.7, respectively. The moment tensor solution of F-net and the spatial distribution of aftershocks indicate that this event was an east-dipping reverse fault type ( Fig. 1a and  b). The earthquake resulted in strong ground motion over Hokkaido, with a maximum seismic intensity of 7 on the JMA scale and a maximum peak ground acceleration of approximately 1800 cm/s/s. The earthquake caused severe damage by strong ground motion and coseismic landslides. Overall, more than 40 people were killed, more than 700 people were injured, and more than 2000 houses were partially or completely destroyed (FDMA 2019).
The center of the island of Hokkaido, where the source area of this earthquake is located, features the Hidaka collision zone, an arc-arc type collision zone between the northeast Japan arc and the Kuril forearc (e.g., Kimura 1996;Kita et al. 2010Kita et al. , 2012. Kita et al. (2012) revealed that beneath the Hidaka collision zone, a low-velocity zone having seismic velocities of crust materials exists in the mantle wedge of the continental plate (the North American Plate) and makes direct contact with the upper surface of the subducting oceanic plate (the Pacific Plate). In the crust beneath the Hidaka collision zone, many earthquakes have occurred, including two recent large Open Access *Correspondence: hkubo@bosai.go.jp 1 National Research Institute for Earth Science and Disaster Resilience, 3-1 Tennodai, Tsukuba, Ibaraki 305-0006, Japan Full list of author information is available at the end of the article reverse-fault-type earthquakes: the 1970 Hidaka earthquake (M JMA 6.7) and the 1982 Urakawa-oki earthquake (M JMA 7.1). The crustal earthquakes within the continental plate have occurred at greater depths compared to ordinary ones in land in the Japanese Islands (HERP 2019). The 2018 Iburi earthquake with a hypocentral depth of approximately 35 km and the subsequent aftershocks also occurred within the crust of the continental plate.
In the conventional framework of strong-motion prediction in Japan, limited attention has been paid to relatively deep crustal earthquakes within the continental plate like the 2018 Iburi earthquake. Therefore, to strengthen the strong-motion prediction methodology, it is important to understand how the rupture of the 2018 Iburi earthquake initiated, grew, and terminated, as well as how the ground motion produced by this earthquake can be simulated using available source and velocity structure models. In this study, we approach these issues using ground-motion records observed by the nationwide strong-motion networks, K-NET and KiKnet, deployed and operated by NIED (e.g., Aoi et al. 2011;NIED 2019b). We first estimate the source process of the event from strong motion waveforms using a multiple-time-window kinematic waveform inversion. Then, using the obtained source model together with a threedimensional (3-D) velocity structure model, we conduct a forward simulation of long-period (< 0.5 Hz) ground motion and compare the result of the simulation with the observations, especially in the Ishikari Lowland. During the 2003 Tokachi-oki earthquake (M JMA 8.0), significant long-period ground motion was observed in the Ishikari Lowland, and the ground-motion characteristics were investigated using a forward waveform simulation (e.g., Hatayama et al. 2007;Aoi et al. 2008).

Data, fault model, and method
The source process is estimated using a multiple-timewindow linear kinematic waveform inversion methodology explained in detail by Sekiguchi et al. (2000) and Suzuki et al. (2010), who followed the approaches of Olson and Apsel (1982) and Hartzell and Heaton (1983).
We use strong-motion waveforms at 20 stations ( Fig. 1a): 14 KiK-net stations in borehole, 4 K-NET stations on ground surface, and 2 F-net (Okada et al. 2004;NIED 2019a) stations in vault. At the F-net stations, we use data observed by velocity-type strong-motion seismographs. Although three components of the seismogram are used at most of the stations, only the horizontal components are used at IBUH03 because the data of the up-down seismograph in the borehole had some problems at the time of the mainshock. Except for data from the F-net stations, where original data are velocity waveforms, the original acceleration waveforms are numerically integrated in the time domain into velocity. The velocity waveforms are band-pass filtered between 0.05 and 0.5 Hz, resampled to 5 Hz, and windowed from 1 s before S-wave arrival for 25 s.
The spatial distribution of the aftershocks shows that the strike angle of the earthquake gradually changed along the strike direction (Fig. 1b). This indicates that it is unreasonable to express the source fault of this event by a single rectangular plane. Consequently, we use a more realistic fault model in the source inversion. Based on the spatial distribution of the aftershocks and the moment tensor solution of F-net, we develop a curved fault model that is approximated by multiple rectangular subfaults ( Fig. 1b and c). The curved fault model is described by multiple planes, each with a width of 20 km, a top depth of approximately 22 km, and a dip angle of 65°, and various strike angles. The dip angle is based on the F-net moment tensor solution. The assumed strike angle of the planes changes gradually along the strike direction: 15° in the northern part, − 20° in the central part, and 40° in the southern part (Fig. 1c). The rupture starting point is set at 42.6908° N, 142.0067° E, at a depth of 37.04 km, which is the hypocenter determined by JMA. The fault model is divided into 15 rectangle subfaults along the strike direction and 10 subfaults along the dip direction.   Kubo et al. Earth, Planets and Space (2020) 72:20 Lowland extending in an NS direction (Fig. 1a); however, the assumed dip angle of this active fault is low at depths greater than 3 km (HERP (Headquarters for Earthquake Research Promotion) 2010), which is inconsistent with the high dip-angle of this earthquake (65° in F-net).
The slip time history of each subfault is discretized using 8 smoothed ramp functions (time windows) progressively delayed by 0.4 s, and having a duration of 0.8 s each. The number of time windows is determined using a trial-and-error process so that the source time function at each subfault can be explained properly. The first time window starting time is defined as the time prescribed by a circular rupture propagation of constant speed, V ftw . Hence, the rupture process and seismic waveforms are linearly related via Green's functions. The slip within each time window at each subfault is derived by minimizing the difference between the observed and synthetic waveforms using a least-squares method. To stabilize the inversion, the slip angle is allowed to vary within ± 45° centered at 107°, which is the rake angle of the F-net moment tensor solution, using the non-negative leastsquares scheme (Lawson and Hanson 1974). In addition, we impose a spatiotemporal smoothing constraint on slips (Sekiguchi et al. 2000). The weight of the smoothing constraint is determined based on Akaike's Bayesian Information Criterion (Akaike 1980).
Green's functions are calculated using the discrete wavenumber method (Bouchon 1981) and the reflection/ transmission matrix method (Kennett and Kerry 1979) with a one-dimensional (1-D) layered velocity structure model. The 1-D velocity structure models are obtained for each station from the 3-D velocity structure model (Fujiwara 2009). Logging data are also referred to for the KiK-net stations. Although a 1-D structure model used for calculating Green's functions differs among stations, the models have a common structure at depth deeper than 13 km where the assumed fault model is located. To consider the rupture propagation effect in each subfault, 25 point sources are uniformly distributed over each subfault to calculate the Green's functions (e.g., Wald et al. 1991). For stations near the fault (IBUH02, HDKH01, HDKH04, IBUH03, and IBUH01), weights that are two times larger than those for the other stations are used. Figure 2a and b shows the total slip distribution of the 2018 Iburi earthquake by map projection and the rupture progression, whereas Additional file 1: Figs. S1 and S2 show the total slip distribution by planar projection and the slip velocity time function of each subfault. The seismic moment is 2.5 × 10 19 Nm (M w 6.9), which is larger than those estimated by F-net and GCMT (1.0 × 10 19 Nm and 1.2 × 10 19 Nm). One likely reason for this is that a portion of slips in the source model are derived to reproduce some part of later phases which are produced by the lateral heterogeneity of the actual velocity structure and which cannot be fully reproduced with the 1-D velocity structure model. One way to approach this issue is further calibration of the velocity structure model for Green's functions. The maximum final slip is 3.8 m. V ftw is set to 1.4 km/s, providing the smallest-misfit solution among all solutions for V ftw varying from 0.8 to 4.0 km/s. The large slips are mostly found at 25-30 km depth in the up-dip direction from the hypocenter. The horizontal extent of this large slip area along the strike direction is approximately 10 km. During the first 6 s, the rupture slowly grows around the hypocenter with small slips. Subsequently, the rupture develops with a large moment release between 6 s and 12 s in the large slip area. The total rupture duration is approximately 20 s.

Results and discussion
The slow rupture growth during the initial moment of this earthquake is supported by the observed waveforms. Figure 2c shows the waveforms recorded at HDKH04 and IBUH01 during the mainshock and of an M JMA 5.4 aftershock that occurred at 6:11 on September 6, 2018 (JST). The latter occurred near the mainshock hypocenter (Fig. 1a), and its focal depth (37.67 km) is practically the same as that of the mainshock. This waveform comparison indicates that although the shape of the main S-wave of the mainshock is similar to that of the aftershock, the main phases of the mainshock arrived several seconds later than those of the aftershock. This suggests that the radiation of seismic waveforms was not significant during the initial moment of this earthquake, which is consistent with the source model derived in this study. Figure 2a shows the spatial distribution of the aftershocks as well as the total slip distribution. Most aftershocks occurred at a depth of approximately 30-35 km, whereas only a few aftershocks occurred within the large slip area. This pattern has also been noted in other crustal earthquakes such as the 1992 Landers earthquake (e.g., Cohee and Beroza 1994), the 1995 south Hyogo (Kobe) earthquake (e.g., Ide et al. 1996), the 1999 Chi-Chi earthquake (e.g., Wu et al. 2001), the 2002 Denali earthquake (e.g., Asano et al. 2005), the 2016 Kumamoto earthquakes (e.g., Kubo et al. 2016), and the 2016 Central Tottori earthquake (e.g., Kubo et al. 2017). Figure 3 shows a comparison between the observed and synthetic waveforms. The waveform fit is satisfactory.

Method and velocity structure model
To evaluate how the ground motion of the earthquake can be reproduced by the current available source and velocity models, we carry out a forward ground-motion simulation at frequencies < 0.5 Hz. We use the 3-D The slip contour interval is 0.8 m. c Comparison of observed waveforms in the mainshock (red) and the M JMA 5.4 aftershock (black) at HDKH04 and IBUH01. These waveforms, which were observed on the ground surface, are low-pass filtered at 0.5 Hz finite-difference method with discontinuous grids (Aoi and Fujiwara 1999;Aoi et al. 2004), in which the grid size is set to 80 m in both the horizontal and vertical directions at depths shallower than 10 km, and three times greater at larger depths.
For the forward simulation, we use the aforementioned source model. For the velocity structure model, J-SHIS v2 3-D velocity model (Fujiwara 2009(Fujiwara , 2012NIED 2019c), which covers the deep subsurface structure from the engineering bedrock (V S = 600 m/s) to the upper crust To validate the performance of the velocity structure model, we first conduct a ground-motion simulation for the M JMA 5.4 aftershock (Fig. 1a), assuming a double-couple point source with the mechanism and the seismic moment referring to F-net, and a smoothedramp type source time function with a 2-s duration. This ground-motion simulation agrees reasonably with the observed ground motion in terms of the velocity waveforms, in the amplitudes of both the direct S-wave pulse and later phases (Additional file 1: Fig.  S4), which implies the appropriateness of the assumed 3-D velocity structure model. The distribution of the peak ground velocity (PGV) is also consistent with the observations, except for the region near stations HKD 129 and HKD180 where the simulation underestimated the observations (Additional file 1: Fig. S5). Note that in this section we use waveforms observed on the ground surface (K-NET and ground-surface seismograph of KiK-net) which are compared with the simulated waveforms at the engineering bedrock (V S = 600 m/s). The effect of site amplification from the engineering bedrock to the ground surface is assumed to be small, as all waveforms are low-pass filtered at 0.5 Hz, although this may not be the case for several stations as discussed later. The PGV in the area around the epicenter is small in contrast to other areas such as the Ishikari Lowland, because this area is located along the node of the S-wave radiation pattern for the aftershock mechanism solution (Fig. 1a), that is, the minimum direction of S-wave radiation. Figure 4 shows the simulated velocity waveforms of the 2018 Iburi earthquake in comparison to the observed waveforms, at selected stations within and around the Ishikari Lowland, along the southern coast of the study area (HDKH06, HKD105, IBUH05, and HKD129), and north of the fault (HKD128, IBUH01, HKD124, HKD182, HKD181, and HKD180) (Fig. 1a). Although the amplitudes of the simulated waveforms are relatively smaller than those of the observed waveforms at most stations, the shape of the S-wave in the simulation is close to that in the observations. The arrival times of the distinctive later phases, which arrive after the S-wave, are consistent with the observations at many stations in the Ishikari Lowland; for instance, at HKD181 at around t = 40-55 s, in the horizontal directions, this consistency is observed (see gray bars in Fig. 4). This suggests that the generation of surface waves in the sedimentary layers could be reproduced by the 3-D wave-propagation simulation. However, at stations HKD129 and HKD180, the amplitudes of the later phases were remarkably underestimated, as they were for the simulation of the M JMA 5.4 aftershock (Additional file 1: Fig. S4), suggesting that the site response of the shallow subsurface structure around these areas, which is not considered in our simulation, may not be negligible even in this frequency range (< 0.5 Hz).

Results and discussion
The spatial distribution of the PGV, or the maximum amplitude of the vector sum of the three-component velocity waveforms, is shown in Fig. 5. The simulated PGV distribution reflects the depth distribution of the seismic bedrock (Additional file 1: Fig. S3), suggesting that the PGV is amplified by the deep subsurface structure. The simulated PGV tends to be consistent with the observations at epicentral distances longer than approximately 40 km, while it tends to be underestimated at shorter distances. Underestimation of the PGV near HKD180 appears both for the mainshock and the M JMA 5.4 aftershock (Additional file 1: Fig. S5), again suggesting the need for considering site responses by the shallow subsurface structure.
The forward waveform simulation of the mainshock, which is presented here, shows that although the overall feature of the observations is reproduced, the amplitude is underestimated. Considering the simulation results of the aftershock in which the observed amplitude was well reproduced at most stations (Additional file 1: Figs. S4 and S5), this result indicates that although the assumed 3-D velocity structure model and the overall feature of the input source model are appropriate, the details in the source model are insufficient. This implies a need for (1) improvement of the source model using the calibration of the velocity structure model for Green's functions, and (2) refinement of the fault model. Previous studies have demonstrated that the uncertainty of input earthquake source largely affects the waveform-simulation results (e.g., Maeda et al. 2016;Iwaki et al. 2017).
In this study, we focused on modeling ground motion at frequencies lower than 0.5 Hz. In future work, it is important to target broadband ground motion to explain the strong ground motion with high seismic intensity (7 as the maximum in JMA scale); this will require precise analysis of the site responses by the shallow subsurface structure. Kubo et al. Earth, Planets and Space (2020) 72:20 Fig. 4 Comparison between observed (black) and simulated waveforms (red) in the forward ground-motion simulation. The number on the left of each waveform denotes the maximum absolute amplitude in cm/s. Distinctive later phases at HKD181 are denoted by gray bars (see text for description)

Conclusions
The 2018 Iburi earthquake was a relatively deep crustal earthquake within the continental plate that caused strong ground motion over a wide area of Hokkaido. Crustal earthquakes at such depth are not accounted for in the current framework of ground-motion prediction in Japan. Therefore, understanding the rupture process and the ground motion caused by the earthquake is significant for improving future predictions. To approach this aim, in this study, we conducted a source inversion and forward ground-motion simulation. The source rupture process of the 2018 Iburi earthquake was first estimated from strong-motion waveforms using a multiple-time-window kinematic waveform inversion. A large slip area with a maximum slip of 3.8 m was found at a depth of 25-30 km in the up-dip direction from the hypocenter. The estimated source model indicated that the rupture slowly grew around the hypocenter with small slips during the first 6 s. Then, the rupture developed with a large moment release between 6 s and 12 s in the large slip area. The large-slip area did not overlap with the high-seismicity area of aftershocks.
Using the obtained source model and the 3-D velocity structure model, J-SHIS v2, we conducted a forward ground-motion simulation in the frequency range of < 0.5 Hz. Although the amplitudes were underestimated at most stations, the arrival times of the distinctive later phases were consistent with the observations at stations of the sedimentary layers of the Ishikari Lowland. The results of the simulation of the mainshock and the aftershock suggest that, although the assumed 3-D velocity structure model and the overall feature of the input source model of the mainshock were appropriate, the details in the source model were insufficient to reproduce the observation of the mainshock.