Source rupture process of the 2016 central Tottori, Japan, earthquake (MJMA 6.6) inferred from strong motion waveforms

The source rupture process of the 2016 central Tottori, Japan, earthquake (MJMA 6.6) was estimated from strong motion waveforms using a multiple-time-window kinematic waveform inversion. A large slip region with a maximum slip of 0.6 m extends from the hypocenter to the shallower part, caused by the first rupture propagating upward 0–3 s after rupture initiation. The contribution of this large slip region to the seismic waves in the frequency band of the waveform inversion is significant at all stations. Another large slip region with smaller slips was found in north-northwest of the hypocenter, caused by the second rupture propagating in the north-northwest direction at 3–5 s. Although the contribution of this slip region is not large, seismic waveforms radiating from it are necessary to explain the later part of the observed waveforms at several stations with different azimuths. The estimated seismic moment of the derived source model is 2.1 × 1018 Nm (Mw 6.1). The high-seismicity area of aftershocks did not overlap with large-slip areas of the mainshock. Two wave packets in the high frequency band observed at near-fault stations are likely to correspond to the two significant ruptures in the estimated source model.Graphical Abstract . .


Background
The 2016 central Tottori earthquake occurred in the central part of Tottori Prefecture, in Chugoku, western Japan, at 14:07 JST on October 21, 2016 (05:07 UTC on October 21, 2016). The Japan Meteorological Agency (JMA) magnitude (M JMA ) is 6.6. The moment magnitude provided by moment tensor analysis of F-net of the National Research Institute for Earth Science and Disaster Resilience (NIED) (Fukuyama et al. 1998) is 6.2. The moment tensor solution of F-net and the spatial distribution of aftershocks determined by Hi-net of NIED indicate that this event was a shallow crustal left-lateral strike-fault type earthquake with a north-northwest (NNW)-southsoutheast (SSE) strike (Fig. 1). This earthquake caused strong ground motions over Tottori and Okayama Prefectures with a maximum seismic intensity of 6-lower on the JMA scale and a maximum peak ground acceleration (PGA) of over 1000 cm/s/s (Fig. 1a).
Tottori Prefecture has been struck by several large crustal earthquakes with severe damage, such as the 1943 Tottori earthquake (M JMA 7.2) (e.g., Kanamori 1972;Kaneda and Okada 2002), the 1983 central Tottori earthquake (M JMA 6.2) (e.g., Nishida 1990;Yoshimura 1994), and the 2000 western Tottori earthquake (M JMA 7.3) (e.g., Fukuyama et al. 2003;Iwata and Sekiguchi 2002). In particular, the 1943 Tottori earthquake caused crippling damage to Tottori City and its surroundings due to strong ground motions. This event killed approximately 1000 people, injured more than 3000 people, and completely or partially destroyed more than 13,000 houses (Usami et al. 2013). Fortunately, the 2016 event did not cause any loss of life, but approximately 30 people were injured and more than 300 houses were completely or partially destroyed (FDMA 2017).
The nationwide strong motion networks, K-NET and KiK-net, deployed and operated by NIED (e.g., Aoi et al. 2011) observed strong ground motions during the 2016 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 central Tottori earthquake. In this study, using the strong motion waveforms, we estimate the source process of the 2016 central Tottori earthquake with a multiple-time-window kinematic waveform inversion. The station coverage for this event was satisfactory ( Fig. 1) and can result in the acquisition of a reliable source model (e.g., Iida et al. 1988;Iida 1990). We compare the coseismic slip distribution of this earthquake with the spatial distribution of aftershocks. We also discuss the relationship between the fault rupture and observed near-fault ground motions. The results of this study can contribute to the accumulation of knowledge on the source and ground-motion characteristics of strikeslip-type crustal earthquakes, and this is useful in improving the ground-motion prediction of future earthquakes.

Data and methods
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 assume a 16 km × 16 km rectangular fault plane model with a strike of 162° and a dip of 88°. The strike and dip angles are based on the F-net moment tensor solution. The rupture starting point is set at 35.3806°N, 133.8545°E, and a depth of 11.58 km, determined by Hi-net. The fault plane is divided into eight subfaults along the strike direction and eight subfaults along the dip direction. The subfault size is 2 km × 2 km. The slip time history of each subfault is discretized using five smoothed ramp functions (time windows) progressively delayed by 0.4 s and having a duration of 0.8 s each. The first time window starting time is defined as the time prescribed by a circular rupture propagation with a 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 −11°, which is the rake angle of the F-net moment tensor solution, using the nonnegative least-squares 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).
We use three-component strong motion waveforms at 15 stations within an epicentral distance of approximately 50 km: seven K-NET stations and eight KiK-net stations (Fig. 1b). K-NET stations are equipped with seismographs on the ground surface, whereas KiK-net stations are equipped with seismographs on the surface and in boreholes. We use the borehole seismograph data at five KiK-net stations, while we use the surface seismograph data at three KiK-net stations (TTRH05, TTRH01, and OKYH10) because their borehole seismograph data had some problems at the time of the 2016 central Tottori earthquake. The original acceleration waveforms are numerically integrated within the time domain into velocity. The velocity waveforms are band-pass filtered between 0.1 and 1.0 Hz with a fourth-order type I Chebyshev filter devised by Saito (1978), resampled to 5 Hz, and windowed from 1 s before S-wave arrival for 10 s.
Green's functions are calculated using the discrete wavenumber method (Bouchon 1981) and the reflection/transmission matrix method (Kennett and Kerry 1979) with a 1-D layered velocity structure model. The 1-D velocity structure model is obtained for each station from the 3-D structure model (Fujiwara et al. 2009). Logging data is also referred to for the KiK-net stations. To consider the rupture propagation effect inside each subfault, 25 point sources are uniformly distributed over each subfault for calculating the Green's functions (e.g., Wald et al. 1991). Figure 2a shows the total slip distribution on the fault. The seismic moment is 2.1 × 10 18 Nm (M w 6.1). V ftw is set to 3.3 km/s, which provides the smallest-misfit solution among all solutions for V ftw varying from 1.6 to 4.0 km/s. A large-slip region, with a maximum slip of 0.6 m, extends from the rupture staring point to the shallower part (R1 in Fig. 2a). Another large-slip region, with a maximum slip of 0.5 m, is located to the NNW of (See figure on previous page.) Fig. 1 a Distribution of peak ground acceleration during the 2016 central Tottori earthquake observed by K-NET and KiK-net. Black star denotes the hypocenter. b Distribution of stations used in this study. Red, blue, and cyan triangles denote K-NET, KiK-net borehole, and KiK-net surface stations, respectively. Black line denotes the assumed fault plane model. Green star denotes the hypocenter of the largest foreshock (October 21, 2016, 12:12 JST, M JMA 4.2). Black circles denote the hypocenters of events (M ≥ 2) from February 1, 2001 to October 20, 2016. Blue circles denote the hypocenters of aftershock within 1 week after the 2016 central Tottori earthquake. These hypocenters were determined by the NIED Hi-net. Indian red, orange, and violet lines denote the source faults of the 1943 Tottori earthquake (M JMA 7.2), the 1983 central Tottori earthquake (M JMA 6.2), and the 2000 western Tottori earthquake (M JMA 7.3), estimated by Kanamori (1972), Nishida (1990), andFukuyama et al. (2003), respectively. Brown, orange, and violet stars denote their hypocenters determined by JMA. Magenta lines denote the surface traces of active faults (AIST 2007). Focal mechanisms represent the F-net moment tensor solutions of the mainshock and the largest foreshock the rupture starting point (R2 in Fig. 2a). The area and the slip values of R1 are larger than those of R2; hence, the released seismic moment on R1 is approximately five times larger than that on R2. Figure 2b, c shows the snapshots of the rupture progression and the slip-velocity time function on each subfault. After the rupture initiation, the rupture mainly propagated from the rupture staring point in the upward direction, and continued for 3 s, causing large slips on R1. Then, the rupture propagated in the NNW direction at 3-5 s, causing large slips on R2. The total rupture duration is approximately 5 s. Figure 3a shows the comparison between the observed and synthetic waveforms. The waveform fit is satisfactory. Figure 3b shows the contribution of the two large regions, R1 and R2 in Fig. 2a, to the synthetic waveforms. At all stations, the contribution of R1 is large. The contribution of R2 is small compared to that of R1. However, Fig. 2 a Total slip distribution on the fault. The slip contour interval is 0.12 m. Vectors denote the direction and the amount of the slip on the hanging wall side. Open star denotes the rupture starting point. Rectangles with blue and magenta broken lines indicate the regions with large slips defined in this study, R1 and R2, respectively. b Rupture progression in terms of slip amount for each 1.0 s interval. The slip contour interval is 0.08 m. c Slip-velocity time function of each subfault. The star denotes the subfault corresponding to the rupture starting subfault. Subfaults emphasized by blue and magenta lines correspond to the subfaults with the largest slip in R1 and R2, respectively. Arrows in the emphasized subfaults denote the peak time of their slip velocities seismic waveforms radiating from R2 are necessary to explain the later part of the observed waveforms at several stations with different azimuths such as TTR007 and OKY006, highlighted with orange bars in Fig. 3b. This result suggests that the slips in R2 are not an artifact. Figure 4 shows the comparison between the slip distribution and the spatial distribution of aftershocks. Many aftershocks occurred around the areas with large slips (>0.4 m) for 1 h after the mainshock, whereas only a few aftershocks occurred within the areas. The aftershocks mainly occurred in the area located approximately 2 km SSE of the mainshock hypocenter (south part of R1) and in the area located approximately 4-7 km NNW of the mainshock hypocenter (north part of R1 and south part of R2). The largest aftershock (October 21, 2016, 14:53 JST, M JMA 5.0) occurred on the south part of R1. In addition, the area of high seismicity within 1 week after the mainshock did not overlap with the large coseismic slip area, although some aftershocks with relatively small magnitude (M < 3) occurred on R1. The feature that aftershocks are not active where coseismic slips are large in the mainshock has also been found in other crustal strike-fault-type earthquakes in Japan: the 1995 south Hyogo (Kobe) earthquake (M JMA 7.2) (e.g., Hirata et al. 1996;Ide et al. 1996), the 2000 western Tottori earthquake (e.g., Ohmi et al. 2002;Semmane et al. 2005), the 2005 west off Fukuoka earthquake (M JMA 7.0) (e.g., Nishimura et al. 2006;Uehira et al. 2006  Next, we focus on the relationship between the source rupture process of the 2016 central Tottori earthquake and ground motions observed at near-fault stations. Velocity waveforms at three near-fault stations, TTR005 (K-NET Kurayoshi), TTRH07 (KiK-net Sekigane), and OKYH10 (KiK-net Kamisaibara) (Fig. 5a), in the low frequency band (0.1-1.0 Hz) are plotted in Fig. 5b. This frequency band is equal to that in the source inversion. The low-frequency waveforms at the near-fault stations have a large pulse at horizontal components. The large pulse is mostly formed by the contribution of R1, and the lowfrequency waveforms from R2 have small amplitudes compared to those from R1 (Fig. 3b). The large amplitude of the low-frequency waveforms from R1 is primarily caused by the large seismic moment release on R1. The difference in the total rupture duration on R1 and R2 could also affect the amplitude difference of the low-frequency waveforms from them; the long rupture duration on R1 could cause a lower dominant frequency of waveforms from R1 than that of waveforms from R2.  Figure 5c shows velocity waveforms at the near-fault stations during the 2016 central Tottori earthquake in the high frequency band (2.0-10 Hz). The high-frequency horizontal waveforms at these stations have two major wave packets, highlighted with orange bars. Figure 5b, c also shows velocity waveforms during the largest foreshock that occurred near the mainshock hypocenter (Figs. 4, 5a) and had a left-lateral strike-fault-type mechanism similar to the mainshock (Fig. 1b). Compared to the mainshock waveforms, the foreshock waveforms are simple: one wave packet in the high frequency band and one pulse in the low frequency band. The arrival time of the first wave packet in the mainshock waveforms is similar to that of the wave packet in the foreshock waveforms.

Results and discussion
The first wave packet in the high-frequency mainshock waveforms corresponds to the large pulse in the lowfrequency waveforms mostly contributed by slips on R1. These indicate that the first wave packet of the high-frequency mainshock waveforms corresponds to the fault rupture on R1. In addition, the foreshock waveforms have no significant later phases, suggesting that the second wave packet of the high-frequency mainshock waveforms should not be attributed to the local underground structure, but to the mainshock fault rupture. Figure 5c shows the observed variation in the arrival-time difference between the two wave packets among the nearfault stations: the arrival-time difference follows the order, TTR005 < TTRH07 < OKYH10. This observation  Fig. 2c). Cyan shade denotes the possible radiation region of the second wave packet. b Observed horizontal velocity waveforms at TTR005 (K-NET Kurayoshi), TTRH07 (KiK-net Sekigane), and OKYH10 (KiK-net Kamisaibara) during the mainshock (black) and the largest foreshock (green) in the low frequency band (0.1-1.0 Hz). The waveform record with borehole observation is shown at TTRH07, while the record with ground surface observation is shown at OKYH10. Broken blue and solid magenta lines denote theoretical arrival times of direct S-waves radiating from the center of the largest-slip subfaults in R1 and R2 (emphasized subfaults in Fig. 2c or colored squares in Fig. 5a) at the peak time of their slip velocities (arrows in the emphasized subfaults in Fig. 2c). c As for Fig. 5b but in the high frequency band (2.0-10 Hz). Orange bars denote two major wave packets in the high-frequency waveforms of the mainshock. d Horizontal envelopes of the high-frequency velocity waveforms. Red line in the envelope at TTR005 denotes its peak time. Light purple zone at TTRH07 denotes the expected range of arrival time of direct S-wave from the radiation region of the second wave packet (cyan shade in Fig. 5a) at the radiation time. Cyan region in the envelope at OKYH10 corresponds to the second wave packet defined in the travel-time analysis suggests that the radiation region of the second wave packet is located to the north of the radiation region of the first wave packet. The potential radiation region of the second wave packet is estimated using the arrival times of the second wave packet in horizontal highfrequency envelope at TTR005 and OKYH10 (Fig. 5d). The details of the analysis are provided in the Additional file 1. The estimated radiation region of the second wave packet (cyan shade in Fig. 5a) is located on the north side of the fault and overlaps R2 in our source model, suggesting that the second wave packet might radiate from R2. We calculate theoretical arrival times of direct S-waves radiating from the center of the largest-slip subfaults in R1 and R2 at the peak time of their slip velocities (Fig. 2c) and plot them in Fig. 5c (blue and magenta lines). Their theoretical arrival times roughly correspond to the observed two wave packets at the near-fault stations and can explain the observed variation in the arrival-time differences among the stations. Thus, it is likely that the two wave packets observed at the near-fault stations correspond to the two significant ruptures at R1 and R2 in our source model. However, these discussions are qualitative and further studies are required to directly estimate the spatiotemporal pattern on radiation strength of high-frequency waveforms.
The waveform appearance in the high frequency band at the near-fault stations differs from that in the low frequency band (Fig. 5): while the low-frequency waveforms have a small amplitude at the time corresponding to the second wave packet, the second wave packet in the high frequency band has a large amplitude compared to the first wave packet. Although the amplitude ratio between the first and second wave packets in the high frequency waveforms differs among stations, the amplitude difference is smaller in the high frequency band than that in the low frequency band. The large amplitude of the second wave packet in high-frequency waveforms is likely due to the dominant frequency of the waveforms from R2 being higher than that of the waveforms from R1 because of the compact fault rupture of R2. In addition, the radiation pattern of body waves, which largely affects the spatial amplitude pattern in the low frequency band, is distorted at frequencies over 1 Hz mainly due to the seismic wave scattering and diffraction within the heterogeneous crust (e.g., Liu and Helmberger 1985;Takemura et al. 2016). Both of these factors lead to the reduction in the amplitude difference shown in the low frequency band; however, they are not enough to explain the larger amplitude of the second wave packet than the first one in the high-frequency waveforms at TTR005. A likely reason is the combination of the non-geometric attenuation with frequency dependence and the short distance from R2. In general, high-frequency waveforms over 1 Hz are likely to be largely affected by the intrinsic and scattering attenuations compared to the low-frequency waveforms. Because the distance from R2 to TTR005 is shorter than that from R1, the amplitude of the high-frequency waveforms radiating from R2 was less affected by the intrinsic and scattering attenuations, causing the large amplitude of the second wave packet in the high-frequency waveforms at TTR005.

Conclusions
In this study, we investigated the source process of the 2016 central Tottori earthquake using a multiple-timewindow kinematic waveform inversion. We used the strong motion waveforms (0.1-1.0 Hz) at 15 stations. In the estimated source model, the seismic moment and the maximum slip are 2.1 × 10 18 Nm (M w 6.1) and 0.6 m, respectively. The source model has two large-slip regions: a region with the maximum slip that extends from the hypocenter to the shallower part (R1), and another with smaller slips located to the NNW of the hypocenter (R2). The released seismic moment on R1 is five times larger than that on R2. Many aftershocks occurred around areas with large coseismic slips, whereas only a few aftershocks occurred within the areas, as observed in other crustal strike-fault type earthquakes in Japan. This earthquake had the foreshock activity near the mainshock hypocenter approximately 2 h before the mainshock.
The contribution of R1 to seismic waves in the low frequency band (0.1-1.0 Hz), which is the analysis frequency band of the waveform inversion, is significant at all stations. Although the contribution of R2 is smaller than that of R1, the seismic waveforms radiating from R2 are necessary to explain the later part of the observed waveforms at several stations with different azimuths. At near-fault stations, two major wave packets were observed in the high frequency band (2.0-10 Hz), and the variation in their arrival times among the near-fault stations indicates that the two wave packets are likely to correspond to the fault ruptures on R1 and R2. The waveform difference in high-and low-frequency bands at the stations can be qualitatively explained using our source model.

Abbreviations
JMA: Japan Meteorological Agency; NIED: National Research Institute for Earth Science and Disaster Resilience; NNW: north-northwest; PGA: peak ground acceleration; SSE: south-southeast.

Additional file
Additional file 1. Estimation of the potential radiation region of the second wave packet observed in the high-frequency waveforms at nearfault stations.