Rupture process of the 2016 Kumamoto earthquake in relation to the thermal structure around Aso volcano

Abstract

We constructed the rupture process model for the 2016 Kumamoto, Japan, earthquake from broadband teleseismic body waveforms (P-waves) by using a novel waveform inversion method that takes into account the uncertainty of Green's function. The estimated source parameters are: seismic moment = 5.1 × 10¹⁹ Nm (Mw = 7.1), fault length = 40 km, and fault width = 15 km. The mainshock rupture mainly propagated northeastward from the epicenter, for about 30 km, along an active strike-slip fault. The rupture propagation of the mainshock decelerated and terminated near the southwest side of the Aso volcano; the aftershock activity was low around the northeastern edge of the major slip area. Our results suggest that the rupture process of the mainshock and the distribution of aftershocks were influenced by the high-temperature area around the magma chamber of Mt. Aso.

Yagi et al. Earth, Planets and Space (2016) 68:118 DOI 10.1186/s40623-016-0492-3

Rupture process ofthe 2016 Kumamoto earthquake inrelation tothe thermal structure aroundAso volcano

Yuji Yagi1*, Ryo Okuwaki2, Bogdan Enescu3, Amato Kasahara2, Ayumu Miyakawa4 and Makoto Otsubo4

Introduction

A large and destructive earthquake (MJMA 7.3) occurred on April 15, 2016 [Coordinated Universal Time (UTC)], in Kumamoto region, Kyushu Island, Japan. The origin time of the earthquake was April 16, 2016, 01:25 [Japan Standard Time (JST)], its epicenter was located at 32.76N, 130.76E, and its hypocentral depth was 12km, as determined by the Japan Meteorological Agency (JMA). Figure 1a shows the largest foreshock (MJMA

6.5) occurred on 14 April, the mainshock, and on- and o-fault aftershock activity within 3days after the 2016 Kumamoto earthquake (epicenters determined by the JMA). The aftershock distribution and a focal mechanism determined in this study (Fig.1a) indicate that the 2016 Kumamoto earthquake occurred along an active strike-slip fault, known as the Futagawa fault (Fig.1a; Research Group for Active Faults of Japan 1991), which belongs to the OitaKumamoto Tectonic Line (OKTL; Fig.1b) (e.g.,

Kamata and Kodama 1994) and is adjacent to the major shear zone in Kyushu (Matsumoto etal. 2015).

The geologic structure, which includes the Aso volcano, locates in the central part of the OKTL and might have contributed to the complex rupture process during fore-shocks, mainshock, and aftershocks. The source model of the mainshock, in particular, may help to understand the seismicity evolution during the 2016 Kumamoto earthquake sequence in its geologic context. Here, we rst determine the robust seismic source process of the 2016 Kumamoto earthquake by applying a novel inversion technique (Yagi and Fukahata 2011a) that takes into account the modeling errors due to the uncertainty of Greens function. Then, we discuss the relationship between the rupture process of the mainshock, seismicity, and the geologic structure around the Aso volcano. Understanding the relationship between seismicity and crustal geophysics (e.g., Hauksson 2011), in particular the crustal thermal structure (Enescu etal. 2009), may help in a better assessment of seismic hazard.

Methods

The spatiotemporal slip-rate distributions of large earthquakes have been estimated by waveform inversion since

*Correspondence: [email protected]

1 Faculty of Life and Environmental Sciences, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8572, JapanFull list of author information is available at the end of the article

2016 The Author(s). This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/

Web End =http://creativecommons.org/licenses/by/4.0/ ), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Yagi et al. Earth, Planets and Space (2016) 68:118

Page 2 of 6

early 1980s (e.g., Hartzell and Heaton 1983). However, the results of inversions conducted for the same earthquake are often dierent from one another (e.g., Beresnev 2003). The modeling errors originated from the accuracy limitations of Greens function are a major problem in seismic source inversion studies and may often bias the inversion result.

To obtain a robust seismic source model for the 2016 Kumamoto earthquake, we applied the inversion formulation developed by Yagi and Fukahata (2011a) to teleseismic P-waves. This inversion formulation takes into account the uncertainty of Greens function and objectively determines the smoothness of the spatiotemporal slip-rate distribution from observed data, using the Akaikes Bayesian information criterion (e.g., Akaike 1980). Such features of the inversion formulation enable us to estimate the complex and irregular rupture process, including the back rupture propagation; the inversion formulation has been applied, for example, to the 2006 Java, Indonesia, tsunami earthquake (Yagi and Fukahata 2011a), the 2011 Tohoku-oki megathrust earthquake

(Yagi and Fukahata 2011b), the 2008 Wenchuan, China, earthquake (Yagi et al. 2012), and the 2015 Illapel, Chile, earthquake (Okuwaki etal. 2016). Yagi and Fukahata (2011a, b) show that if we neglect the data covariance components, the slip-rate distribution is distorted by modeling errors originated from the uncertainty of Greens function.

We have used the vertical component teleseismic P-wave data recorded at 27 broadband seismic stations (Fig. 2a). In general, we can obtain the robust slip-rate distribution of large earthquakes from teleseismic P-wave data, which are insensitive to the assumed structure (Yagi etal. 2004). We removed the seismograph response from the original waveforms to convert them into velocity waveforms and then re-sampled the waveform data at 0.3-s intervals. In order to mitigate aliasing and low frequency noise, a 0.001- to 0.56-Hz Butterworth band-pass lter was applied before re-sampling. The rst arrival of P-phase on each original waveform was manually picked so that it corresponds to the initial rupture break.

Yagi et al. Earth, Planets and Space (2016) 68:118

Page 3 of 6

We assumed the rupture occurring on a single fault plane of (strike, dip) = (234, 64) and the fault area of 58km length and 20km width. The geographical coordinates of the initial break point are the same as for the JMA epicenter. The fault plane was slightly modied from that obtained by moment tensor inversion of teleseismic body waves to be consistent with the geometry and location of the Futagawa fault (Headquarters for Earthquake Research Promotion, http://www.jishin.go.jp/main/index%e2%80%90e.html

Web End =http://www.jishin. http://www.jishin.go.jp/main/index%e2%80%90e.html

Web End =go.jp/main/indexe.html ; last accessed on April 28, 2016). For the same reason, the break point depth (i.e., hypo-central depth) was taken at 9 km. Assumed dip angle was shallower than the estimated value of the regional seismic network (National Research Institute for Earth Science and Disaster Resilience, Japan, http://www.fnet.bosai.go.jp/top.php?LANG=en

Web End =http://www.fnet. http://www.fnet.bosai.go.jp/top.php?LANG=en

Web End =bosai.go.jp/top.php?LANG=en ; last accessed on June 13, 2016). We performed inversions using several dip angles and found that the center of the maximum slip area was stable (Additional le 1: Figure A2). We imposed no specic constraints on the rake angles. As the structure model near the source, we used CRUST 1.0 (Additional le 1: Table A1; Laske et al. 2013). We adopted a slip-rate function represented as linear B-splines with a grid interval of 0.3 s and having a duration of 15 s on each fault patch. The assumed rupture duration was taken as 18 s. The theoretical Greens functions of each source node were calculated with the method of Kikuchi and

Kanamori (1991). The maximum rupture front velocity, which denes the rupture onset time at each spatial node, was set to 2.4 km/s. The center of the maximum slip area was stable even when we have assumed various maximum rupture front velocities (Additional le1: Figure A3). We also examined two kinds of initial slip model: uniform slip model and vertical striped slip models in the same way as Yagi and Fukahata (2011a), and found that exactly the same solution was obtained for both of the initial slip models.

Results anddiscussion

The total slip distribution, the moment-rate function, and the focal mechanism are given in Fig.1a. We found that the rupture mainly propagated toward northeast from the epicenter. The rupture direction is consistent with the characteristics of the mainshock waveforms (Fig.2b), which are shorter toward northeast and broader toward southwest. The area of large slip, with a maximum of 5.7m, is centered about 10km northeast from the epicenter. The eective slip area extends about 40km long and 15km wide. The average rake angle was 148, which is approximately consistent with the fault motion of historical earthquakes estimated by geological surveys (Research Group for Active Faults of Japan 1991). The total seismic moment was 5.1 1019 Nm (Mw = 7.0), which is comparable to the Global Centroid Moment

Yagi et al. Earth, Planets and Space (2016) 68:118

Page 4 of 6

Tensor solution of 4.51019Nm (http://www.globalcmt.org

Web End =http://www.globalcmt. http://www.globalcmt.org

Web End =org ; last accessed on April 28, 2016).

Figure3 shows the distribution of foreshocks, the spatiotemporal slip rate of the mainshock, and the aftershocks. The foreshocks started in the northern part of an active strike-slip fault, known as the Hinagu fault (Fig. 3a), and expanded northeastward and southwest-ward along the Hinagu fault. The rupture of the main-shock started near the junction of the Hinagu and the

Futagawa fault zones. Such a foreshockmainshock rupture pattern is well known for complex fault zones (e.g., King and Nabelek 1985), suggesting that the foreshock activity near the junction between the Hinagu and Futagawa faults contributed to the initiation of the mainshock rupture.

The slip-rate distribution shows that the rupture propagated to southwest during the rst 5s. At 5s after the origin time, the main rupture started to propagate toward

Yagi et al. Earth, Planets and Space (2016) 68:118

Page 5 of 6

northeast. The slip rate started decreasing from 15s and terminated at about 17s after the initial break. The on-fault aftershocks occurred along the Futagawa fault zone, and the o-fault aftershocks were triggered along the OKTL (Figs.1a, 3c).

To better understand the complex relationship between the geologic structure, the rupture process of the main-shock, and the aftershock activity, we compare in Fig.4 the land elevation (Fig.4a), aftershock density (Fig.4b), and shear stress change along the fault strike (Fig.4c). We calculated the shear stress change corresponding to the average slip direction on the assumed fault plane using the formulation for a static dislocation in a homogeneous half-space (Okada 1992). The maximum stress drop was 9.5MPa (Fig.4c) near the largest slip area (Fig.1a). On the other hand, we notice a zone of stress increase (stress rise values up to 5.6 MPa) immediately northeast of the large stress drop zone. The aftershock activity (Fig.4b) is highest in the hypocentral region [Fig.4b: (5 to +10km), along strike], next to the area of large stress drop during the mainshock. Note that the enhanced aftershock activity here is likely the result of superposed aftershocks of the mainshock and largest (MJMA 6.5)

foreshock (Fig. 3a, c). The along-strike aftershock density decreases gradually toward northeast [0 to 20km, along strike]. The moderate aftershock activity taking place at the same along-strike distances [Fig.4b: (5 to 15km), along strike] as the maximum shear stress drop reects mainly deeper aftershock activity (i.e., below the large slip area, see Figs.1a, 4c). Further to the northeast, there is a quite sudden decrease in the aftershock activity toward very low levels [Fig. 4b: (20 to 30 km), along strike], followed by an along-strike increase that reaches a maximum [Fig. 4b: ~(40 km), along strike]

around the northeastern edge of the Aso volcano. Note that it is difficult to discuss the relationship between this northeastern-most increase in seismicity and the shear stress change due to small slip near the northeastern edge (Figs.1a, 4c), since the uncertainty of slip distribution increases with distance from the hypocenter (Yagi and Fukahata 2011a). However, it is clear that the aftershock activity is low around the rupture deceleration area, where the shear stress was increasing due to the coseismic slip of the mainshock, and thus, we would have expected an increased aftershock activity. In general, high-temperature areas promote plastic deformation and inhibit dynamic rupture propagation during earthquakes (e.g., Scholz 1998). A tomography study conducted in the region shows that the low seismic velocity sphere, located shallower than 10km beneath the Aso volcano, relates to the existence of a magma chamber (Sudo and Kong 2001). The thermal structure model determined by magnetic data (Okubo and Shibuya 1993) shows that the

high temperature (over approximately 500C) reached to about 3km beneath the Aso caldera. The high-temperature area around the magma chamber of Mt. Aso may have contributed to the termination of the rupture during the mainshock and the low aftershock activity around the Aso volcano.

Conclusions

We estimated the robust spatiotemporal slip-rate distribution of the 2016 Kumamoto earthquake using the waveform inversion of teleseismic P-wave data and examined the 2016 Kumamoto earthquake sequence, including the foreshock and aftershock activity. The foreshocks started

Yagi et al. Earth, Planets and Space (2016) 68:118

Page 6 of 6

in the northern part of the Hinagu fault, and the rupture of the mainshock initiated at the junction between the Hinagu and the Futagawa fault zones. The dynamic rupture of the mainshock mainly propagated 30 km north-eastward from the epicenter, along the Futagawa fault. The rupture of the mainshock decelerated and terminated near the southwest side of the Aso volcano, and aftershock activity was low around the northeastern edge of the major slip area (where shear stress was increasing due to the coseismic slip of the mainshock). The high-temperature area around the magma chamber of Mt. Aso may contribute to the termination of the rupture during the mainshock and peculiarly low aftershock activity of the 2016 Kumamoto earthquake, around Mt. Aso.

References

Akaike H (1980) Likelihood and the Bayes procedure. In: Bernardo JM, DeGroot MH, Lindley DV, Smith AFM (eds) Bayesian Statistics. University Press, Valencia, pp 143166Beresnev IA (2003) Uncertainties in nite-fault slip inversions: to what extent to believe? (a critical review). Bull Seismol Soc Am 93:24452458. doi:http://dx.doi.org/10.1785/0120020225

Web End =10.1785/0120020225 Enescu B, Hainzl S, Ben-Zion Y (2009) Correlations of seismicity patterns in

Southern California with surface heat ow data. Bull Seismol Soc Am 99(6):31143123. doi:http://dx.doi.org/10.1785/0120080038

Web End =10.1785/0120080038 Farr TG, Rosen PA, Caro E, Crippen R, Duren R, Hensley S, Kobrick M, Paller M,

Rodriguez E, Roth L, Seal D, Shaer S, Shimada J, Umland J, Werner M, Oskin M, Burbank D, Alsdorf D (2007) The shuttle radar topography mission. Rev Geophys 45:RG2004. doi:http://dx.doi.org/10.1029/2005RG000183

Web End =10.1029/2005RG000183 Hartzell S, Heaton T (1983) Inversion of strong ground motion and teleseismic waveform data for the fault rupture history of the 1979 Imperial Valley, California, earthquake. Bull Seismol Soc Am 73:15531583Hauksson E (2011) Crustal geophysics and seismicity in southern California.

Geophys J Int 186:8298. doi:http://dx.doi.org/10.1111/j.1365-246X.2011.05042.x

Web End =10.1111/j.1365-246X.2011.05042.x Kamata H, Kodama K (1994) Tectonics of an arc-arc junction: an example from

Kyushu Island at the junction of the Southwest Japan Arc and the Ryukyu Arc. Tectonophysics 233:6981Kikuchi M, Kanamori H (1991) Inversion of complex body wavesIII. Bull

Seismol Soc Am 81:23352350King G, Nabelek J (1985) The role of fault bends in the initiation and termination of earthquake rupture. Science 228:984987Komazawa M (2013) Gravity grid data of Japan. In: Survey Geological, Geological Survey of Japan (eds) Gravity database of Japan DVD edition Digital Geoscience Map P-2. Geological Survey of Japan, AIST, TsukubaLaske G, Masters G, Ma Z, Pasyanos M (2013) Update on CRUST1.0A 1-degree

Global Model of Earths Crust. Geophys Res Abstracts 15:EGU2013-2658 Matsumoto S, Nakao S, Ohkura T, Miyazaki M, Shimizu H, Abe Y (2015) Spatial heterogeneities in tectonic stress in Kyushu. Earth Planets Space, Japan and their relation to a major shear zone. doi:http://dx.doi.org/10.1186/s40623-015-0342-8

Web End =10.1186/s40623-015-0342-8 Okada Y (1992) Internal deformation due to shear and tensile faults in a half-space. Bull Seism Soc Am 82:10181040Okubo Y, Shibuya A (1993) Thermal and crustal structure of the Aso volcano and surrounding regions constrained by gravity and magnetic data, Japan. J Volcanol Geotherm Res 55:337350Okuwaki R, Yagi Y, Arnguiz R, Gonzlez J, Gonzlez G (2016) Rupture process during the 2015 Illapel, Chile earthquake: zigzag-along-dip rupture episodes. Pure appl Geophys 173:10111020. doi:http://dx.doi.org/10.1007/s00024-016-1271-6

Web End =10.1007/ http://dx.doi.org/10.1007/s00024-016-1271-6

Web End =s00024-016-1271-6 Research group for active faults of Japan (1991) Active faults in Japan. Sheet maps and inventories, Revised Edition. University of Tokyo Press, Tokyo Scholz CH (1998) Earthquake and friction laws. Nature 391:3742. doi:http://dx.doi.org/10.1038/34097

Web End =10.1038/34097 Sudo Y, Kong L (2001) Three-dimensional seismic velocity structure beneath

Aso Volcano, Kyushu, Japan. Bull Volcanol 63:326344. doi:http://dx.doi.org/10.1007/s004450100145

Web End =10.1007/ http://dx.doi.org/10.1007/s004450100145

Web End =s004450100145 Wessel P, Smith WHF (1998) New, improved version of generic mapping tools released. EOS Trans Am Geophys Union 79:579. doi:http://dx.doi.org/10.1029/98EO00426

Web End =10.1029/98EO00426 Yagi Y, Fukahata Y (2011a) Introduction of uncertainty of Greens functioninto waveform inversion for seismic source processes. Geophys J Int 186:711720. doi:http://dx.doi.org/10.1111/j.1365-246X.2011.05043.x

Web End =10.1111/j.1365-246X.2011.05043.x Yagi Y, Fukahata Y (2011b) Rupture process of the 2011 Tohoku-oki earthquake and absolute elastic strain release. Geophys Res Lett 38:16. doi:http://dx.doi.org/10.1029/2011GL048701

Web End =10.1029/ http://dx.doi.org/10.1029/2011GL048701

Web End =2011GL048701 Yagi Y, Mikumo T, Pacheco J, Reyes G (2004) Source rupture process of the

Tecoman, Colima, Mexico earthquake of January 22, 2003, determined by joint inversion of teleseismic body wave and near source data. Bull Seism Soc Am 94:17951807Yagi Y, Naoki N, Kasahara A (2012) Source process of the 12 May 2008

Wenchuan, China, earthquake determined by waveform inversion of teleseismic body waves with a data covariance matrix. Earth Planets Space 64:e13e16. doi:http://dx.doi.org/10.5047/eps.2012.05.006

Web End =10.5047/eps.2012.05.006

Additional le

Abbreviations

JMA: Japan Meteorological Agency; JST: Japan Standard Time; OKTL: Oita Kumamoto Tectonic Line; SRTM: Shuttle Radar Topography Mission; UTC: Coordinated Universal Time.

Authors contributions

YY led and designed the whole research and drafted the manuscript. RO carried out the inversion analysis, generated the gures and Additional le, and contributed to the examination of the results. BE and AK contributed to the discussion of the results. AM contributed to generate the contents of Fig. 1b. AM and MO contributed to the interpretation of the results related to the local and regional geology. All authors discussed the results and commented on the manuscript. All authors read and approved the nal manuscript.

Author details

1 Faculty of Life and Environmental Sciences, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8572, Japan. 2 Graduate School of Lifeand Environmental Sciences, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8572, Japan. 3 Department of Geophysics, Faculty of Science, Kyoto University, Kita-Shirakawa Oiwake-cho, Sakyo-ku, Kyoto 606-8502, Japan.

4 Geological Survey of Japan, AIST, 1-1-1 Higashi, Tsukuba, Ibaraki 305-8568, Japan.

Acknowledgements

Waveform data from Berkeley Digital Seismic Network, Global Seismograph Network, and IRIS/Singapore, Singapore National Network, were accessed through the Incorporated Research Institutions for SeismologyData Management Center. Figures were generated with the Generic Mapping Tools (Wessel and Smith 1998). We thank Yasukuni Okubo and Yujiro Ogawa for useful information. We also thank Chen Ji, Sebastiano DAmico, and the Editor Haruo Horikawa and the Chief Editor Yasuo Ogawa for helpful comments and advices. Authors acknowledge support from the Japan Society for the Promotion of Science (JSPS) KAKENHI Grant 16K05529 (to YY) and 26240004 (to BE).

Competing interests

The authors declare that they have no competing interests.

Received: 5 May 2016 Accepted: 15 June 2016