In this study, the rupture process of the main shock of the Kumamoto earthquake, particularly the generation of strong ground motions in the frequency range relevant to structural damage, was investigated based on the inversion of strong ground motions. Strong-motion records in the near-source region were mainly utilized because the authors were interested in the generation mechanism of damaging ground motions in the near-source region. Empirical Green’s functions (EGFs) were applied to avoid uncertainty in the subsurface structure model. Four cases of inversions with different combinations of small events were used to investigate the dependence of the inversion results on the selection of the small events. It was found that the dependence of the final slip distribution and peak slip velocity distribution on the selection of the EGF events is small. The results clearly indicate that a region of significantly large slip and slip velocity existed approximately 15 km northeast of the hypocenter. However, no “asperity” was observed between the hypocenter and Mashiki. Thus, it is not appropriate to conclude that the large-amplitude pulse-like ground motion in Mashiki was generated by the forward-directivity effect associated with the rupture of an asperity. As far as the source effect is concerned, the ground motion in Mashiki cannot be interpreted as the worst case scenario. On the other hand, the rupture of the “asperity” 15 km northeast of the hypocenter should have caused significantly large ground motions in regions close to the asperity. The significant damage of highway bridges in the region can potentially be attributed to the rupture of the asperity. The result of this study was compared with an inversion result obtained from numerical Green’s functions for a layered half-space. The two results share the main features in spite of the difference of the Green’s functions and stations used. Therefore, it can be concluded that these two source models capture the main features of the rupture process of the earthquake.

Graphical abstract

In this study, the rupture process of the main shock of the Kumamoto earthquake, particularly the generation of strong ground motions in the frequency range relevant to structural damage, was investigated based on the inversion of strong ground motions. Strong-motion records in the near-source region were mainly utilized because the authors were interested in the generation mechanism of damaging ground motions in the near-source region. Empirical Green’s functions (EGFs) were applied to avoid uncertainty in the subsurface structure model. Four cases of inversions with different combinations of small events were used to investigate the dependence of the inversion results on the selection of the small events. It was found that the dependence of the final slip distribution and peak slip velocity distribution on the selection of the EGF events is small. The results clearly indicate that a region of significantly large slip and slip velocity existed approximately 15 km northeast of the hypocenter. However, no “asperity” was observed between the hypocenter and Mashiki. Thus, it is not appropriate to conclude that the large-amplitude pulse-like ground motion in Mashiki was generated by the forward-directivity effect associated with the rupture of an asperity. As far as the source effect is concerned, the ground motion in Mashiki cannot be interpreted as the worst case scenario. On the other hand, the rupture of the “asperity” 15 km northeast of the hypocenter should have caused significantly large ground motions in regions close to the asperity. The significant damage of highway bridges in the region can potentially be attributed to the rupture of the asperity. The result of this study was compared with an inversion result obtained from numerical Green’s functions for a layered half-space. The two results share the main features in spite of the difference of the Green’s functions and stations used. Therefore, it can be concluded that these two source models capture the main features of the rupture process of the earthquake.

A series of damaging earthquakes hit the Kumamoto and Oita prefectures, Kyushu, Japan, including the Mw6.1 foreshock (April 14, 21:26 JST) and the Mw7.1 main shock (April 16, 1:25 JST). As of July 12, it involved 1879 perceptible earthquakes (Japan Meteorological Agency 2016a). The entire sequence was named “the 2016 Kumamoto earthquake” by the Japan Meteorological Agency (JMA; Japan Meteorological Agency 2016b). Eighty-one people were killed because of the earthquakes according to the Fire and Disaster Management Agency, Japan (as of July 29).

The Mw7.1 main shock was recorded not only by permanent strong-motion networks, such as K-NET (Kinoshita 1998; Aoi et al. 2004), KiK-net (Aoi et al. 2000, 2004), and the JMA network (Nishimae 2004), but also by a temporary seismic array that was deployed after the Mw6.1 foreshock in Mashiki (see Fig. 1), Kumamoto Prefecture, located in the near-source region (Hata et al. 2016). The importance of their observation is that they successfully recorded strong ground motions during the main shock in a heavily damaged zone in which the ratio of totally collapsed wooden houses exceeded 50% (Kikuchi and Tanaka 2016), although the damage could have partly been due to the Mw6.1 foreshock. The ground motions observed during the main shock involved a large-amplitude pulse with an approximate period of 1 s and a peak velocity of 1.59–1.79 m/s (Hata et al. 2016). In terms of spectral acceleration, the observed ground motions exceeded the largest observed ground motions of the 1995 Kobe earthquake in the period around 1 s. It is important to understand what kind of rupture process contributed to the generation of such strong ground motions.

Fig. 1

Fault plane assumed in the inversion and locations of strong-motion stations plotted on a topological map by the Geospatial Information Authority of Japan. The long rectangle oriented in N52E direction indicates the fault plane assumed in the inversion. The star is the rupture starting point assumed in the inversion, which is the hypocenter. The open rectangles indicate the hypocenters of the small events used in the analysis (EGF1–EGF4). The solid triangles are the locations of strong-motion stations used in the inversion. The open triangles represent the locations of strong-motion stations used in “Discussion” section

In addition to the damage to buildings, damage to highway bridges was also significant in the near-source region of the Kumamoto earthquake. Although some bridges were clearly affected by landslides, several bridges, including the Ogiribata Bridge (see Fig. 1), could have been simply damaged by strong ground motions (Nikkei Construction 2016). To understand the true damage mechanism of the bridges, it is necessary to estimate the strong ground motions at the bridge sites. For this purpose, again, it is necessary to understand the rupture process of this earthquake especially when the bridge site is located close to the fault, because the rupture process of the earthquake easily affects the strong ground motions near the fault.

Therefore, the rupture process of the main shock of the Kumamoto earthquake was investigated in this study based on the inversion of strong ground motions, in particular the generation of strong ground motions in the frequency range relevant to structural damage. Strong-motion records in the near-source region were mainly used in the study because the authors were interested in the generation mechanism of damaging ground motions in the near-source region.

Empirical Green’s functions (EGFs) were used in the waveform inversion in this study. As discussed in Nozu (2007), the main advantage of the EGFs over numerical Green’s functions is that they allow us to avoid uncertainty in the subsurface structure model used in calculating the numerical Green’s functions, which could affect the estimated rupture process; even complex 2D or 3D effects of the ground motions can be naturally incorporated using EGFs. The 3D nature of the subsurface structure is particularly significant in the near-source region of the Kumamoto earthquake, as represented by the existence of the Aso caldera (Fig. 1). In fact, the effect of the 3D subsurface structure is quite noticeable in the observed waveforms, as discussed later in the article. Therefore, the use of EGFs could be especially advantageous for this particular earthquake.

One of the interesting questions associated with the inversion results obtained using EGFs is to what extent the inversion results are dependent on the selection of the EGFs. To answer this question, four cases of inversions with different combinations of EGFs were used in this study. The results were then compared in terms of final slip distributions and peak slip velocities.

Because the purpose of the study was to reveal the rupture process of the earthquake that caused damaging ground motions in the near-source region, near-source stations were used in the analysis. In particular, five K-NET stations (KMM004, KMM005, KMM006, KMM009, KMM011) and three KiK-net stations (KMMH06, KMMH14, KMMH16) were selected (solid triangles in Fig. 1). The KiK-net station KMMH16 is located in Mashiki. Although it is located outside of the heavily damaged zone (Hata et al. 2016), it is still within 1 km from the temporary stations deployed in the heavily damaged zone. The observed velocity waveforms at KMMH16, both at the surface and in the borehole, were also characterized by a large-amplitude pulse with an approximate period of 1 s, similar to those observed at the temporary stations, although the peak ground velocity at the surface was 20–25% smaller than at the temporary stations. Therefore, information about the effect of the earthquake rupture on the damaging ground motions in Mashiki could be obtained from including KMMH16 in the inversion analysis.

Although the use of near-source stations is advantageous in constraining the detailed spatiotemporal distribution of slip on the fault plane, its downside is that the inversion results could be affected by soil nonlinearity, because the amplitude of the strong ground motions tends to be large in the near-source region. To minimize the effect of soil nonlinearity, borehole records, rather than surface records, were used at the KiK-net stations. However, surface records were used at KMMH06 because the S/N ratio of the borehole records was not sufficient for small events. At the K-NET stations, surface records were used. The potential effects of soil nonlinearity on the results will be discussed in “Discussion” section.

Figure 2 shows the EW and NS components of the main shock velocity waveforms (0.2–2 Hz) at the target stations (gray and black traces), which were used in the inversion analysis. The waveforms were obtained by integrating the accelerograms in the frequency domain. The target frequency range was determined so that it covers the 1-Hz energy predominant at KMMH16. Time sections with a duration of 15 s, including the S waves, were used in the waveform inversion (black traces in Fig. 2). The waveforms at KMM004 and KMMH06, which are both located in the Aso caldera, are characterized by the predominance of the later phases with long duration. These observations indicate complex 3D effects of the subsurface structure at these sites. The use of empirical Green’s functions is therefore advantageous for the present analysis.

Fig. 2

Comparison of the observed and synthetic velocity waveforms (0.2–2 Hz) for Case 1. The gray and black traces show the observed waveforms. The black traces correspond to the 15-s time sections used in the waveform inversion. The red traces show the synthetic waveforms

In the empirical Green’s function method, it is important to use small events which share the path and site effects with the large event. Because the main shock fault plane of the 2016 Kumamoto earthquake is relatively large and the target stations are located relatively close to the fault plane (Fig. 1), it might not be appropriate to represent the contributions of the entire fault plane with the records of only one small event. Thus, the main shock fault plane was divided into two parts: the 16-km-long western part and the 24-km-long eastern part. The different small events were then assigned to those parts as were done in previous studies (e.g., Nozu 2007). It should be noted that the station KMMH16 approximately corresponds to the boundary. Four cases of inversions (Cases 1–4) with different combinations of small events were used to investigate the dependence of the inversion results on the selection of the small events (EGF1–EGF4; Table 1). The parameters of those events are listed in Table 2. In Case 1, the events at 0:50 JST (EGF1) and 15:27 JST (EGF2) on April 15 were selected for the western and eastern parts of the fault plane, respectively, because these events met the following conditions:

1.

The hypocenters are located close to the respective parts of the main shock fault plane (Fig. 1).

2.

The mechanisms estimated by F-net (Fukuyama et al. 1996) are similar to that of the main shock; they are right-lateral earthquakes occurring on a near-vertical fault (Table 2).

3.

MJ ≥ 4.0. This condition was imposed because the records of smaller earthquakes tend to exhibit insufficient S/N ratios at low frequencies.

4.

The records are obtained at all target stations.

5.

The Fourier phase characteristics of the small event resemble those of the main shock (Nozu 2007; Nozu and Irikura 2008).

The combination selected for Case 1 is exactly the same as in the preliminary analysis (Nozu 2016). Case 1 is regarded as the reference case. Subsequently, alternative small events were searched. It was impossible to find a small event that satisfies all of the conditions 1–5 in the western part of the fault. Therefore, the event on April 16 (4:51 JST; EGF3), which satisfies conditions 1, 3, 4, and 5, was selected. Although the moment tensor solution of the F-net was not available for the earthquake, the hypocenter location was appropriate (Fig. 1) and the Fourier phase characteristics were preferable compared to that of other earthquakes, especially at KMM006. It was also impossible to find a small event that satisfies all of the conditions 1–5 in the eastern part of the fault. Therefore, the event at 7:29 JST on April 15 (EGF4), which satisfies the conditions 1, 3, 4, and 5, was selected. The moment tensor solution of the F-net was not available for the earthquake. The EW and NS components of the velocity waveforms of the EGF events (0.2–2 Hz) at the target stations were obtained by integrating the accelerograms in the frequency domain in the same way as the main shock waveforms; they were used as empirical Green’s functions.

Linear least-squares waveform inversion (Nozu 2007; Nozu and Irikura 2008) was adopted. The inversion follows the multi-time-window approach, which was proposed by Hartzell and Heaton (1983), although empirical Green’s functions were used in the present analysis. The method can be summarized as follows:

First, the fault plane of the large event is divided into NL × NW subfaults. At each subfault, the moment rate function after passage of the first-time-window triggering front is then expressed as a convolution of an impulse train, with the moment rate function of the small event allocated to that particular subfault as follows:

where Mij(t) is the moment rate function at the ij subfault (t = 0 corresponds to the passage of the first-time-window triggering front), m(t) is the moment rate function of the small event (t = 0 corresponds to the rupture time of the small event), ND is the discretization number along the time axis, Δτ is the interval of discretization, and wijk is the ratio of the moment release at the ij subfault in the time interval [(k − 1)Δτ, kΔτ] with respect to the moment of the small event, which is the unknown parameter of the inversion. At first glance, it may seem that Mij(t) in Eq. (1) involves a “gap” when the source duration of the small event is shorter than the time interval. However, it should be noted that the present study targets the low-frequency components (<2 Hz). Therefore, only the low-pass-filtered version of Eq. (1) is used in this study. As an example, when the time interval is 0.25 s, as it is the case in this study, the low-pass-filtered (at 2 Hz) version of the impulse train in Eq. (1) does not involve any gap. Thus, the source duration of the small event can be shorter than the time interval. On the other hand, if the source duration of the small event is much greater than the time interval, the resolution of the inversion is degraded. It would be difficult to estimate the moment release for each time window. However, the corner frequencies of the EGF events are sufficiently high in the present study, which will be discussed later. Therefore, we can estimate the moment release for each time window.

where V(t) is the velocity waveform for the large event (t = 0 corresponds to the beginning of the record), v(t) is the velocity waveform for the small event (t = 0 corresponds to the beginning of the record), ra is the hypocentral distance of the small event, rij is the distance from the ij subfault to the observation station, VS is the shear wave velocity in the source region, ξij is the distance from the rupture starting point of the large event to the ij subfault, VR is the velocity at which the first-time-window triggering front propagates, ta is the record time of the small event, ta0 is the rupture time of the small event, tm is the record time of the large event, and tm0 is the rupture time of the large event. In this formulation, it is assumed that a subfault and the small event associated with the subfault share the same rake angle. Once the unknown parameters are determined and combined with the moment of the small event, the spatiotemporal distribution of the moment release can be obtained. The spatiotemporal distribution of the slip can then be obtained by dividing it with the rigidity. As shown in Eq. (4), it is not necessary to detect the initial S phase in the records of the large event in this formulation. Therefore, the method is advantageous for earthquakes with an ambiguous initial S phase.

The relation between the size of the small event, size of the subfault, and frequency component to be analyzed is as follows: In this formulation, a sufficiently small event that can be regarded as a point source is used (condition 1). Similarly, a sufficiently small subfault, which can be regarded as a point source, is utilized (condition 2). Based on condition 1, sufficiently low-frequency components should be used in the analysis in light of the size of the small event. In general, the corner frequency involved in the source spectrum of an earthquake is related to the finiteness of the fault size and the rise time (e.g., Aki and Richards 2002); the effect of the fault finiteness is recognizable at frequencies higher than the corner frequency. Therefore, condition 1 requires the use of frequency components lower than the corner frequency of the small event. Condition 2 is the same as that for the waveform inversion with numerical Green’s functions. It requires the use of sufficiently small subfaults, depending on the frequency components to be analyzed. Thus, conditions 1 and 2 require that both the area of the small event and that of the subfault are sufficiently small. However, it is not necessary that those areas coincide with each other. In this respect, the present method differs from the empirical Green’s function method of Irikura (1986).

Non-negative least-squares solutions are obtained using the algorithm of Lawson and Hanson (1974). In addition, constraints are imposed to minimize the second-order derivative of the slip, with respect to both time and space (Nozu 2007). The constraints are imposed as follows:

In terms of the second-order derivative with respect to space, it should be noted that Eq. (5) is also applied at the outer boundaries of the fault (i.e., for i = 1, i = NL, j = 1, and j = NW). At the outer boundaries, Eq. (5) is used with the assumption that w = 0 outside of the fault plane. This constraint is imposed with the intention of obtaining a slip distribution that converges to zero at the outer boundaries. At the boundaries of regions associated with different EGF events, the absolute moment at each subfault (rather than the moment ratio to the EGF event) is assumed to be continuous; wijk is replaced by wijkmij in Eq. (5), where mij is the moment of the EGF event associated with the ij subfault.

In terms of the second-order derivative with respect to time, Eq. (6) is applied for k = ND but not for k = 1, with the intention of allowing high slip velocity just after passage of the first-time-window triggering front.

For the application to the main shock of the 2016 Kumamoto earthquake, a fault plane that includes the JMA hypocenter with a length of 40 km and width of 20 km was assumed (Fig. 1). The strike and dip angles were initially set to 226° and 84°, respectively, referring to the F-net moment tensor solution (Table 2). The strike angle was then changed to 232° to be more consistent with the surface fault trace.

The fault was divided into 20 × 10 subfaults. It was assumed that the first-time-window triggering front starts at the JMA hypocenter and propagates radially at a constant velocity. The velocity will be referred to as the “first-time-window triggering velocity” in this article. The choice of the velocity will be discussed later. In Eq. (1), an impulse train that spans 3.0 s and consists of 12 impulses at equal time intervals of 0.25 s was used. Thus, the height of each impulse corresponds to the ratio of the moment release during 0.25 s, with respect to the moment of the EGF event. The ratio was determined using the inversion. In terms of the appropriateness of the total time window of 3.0 s, the authors analyzed another case of inversion in which the total time window was extended to 4.0 s. As a result, it was confirmed that the final slip distribution was not significantly affected by this modification.

The shear wave velocity in the source region was assumed to be 3.55 km/s according to Fukuyama et al. (1998).

The moments of EGF1 and EGF2 determined by the F-net (Fukuyama et al. 1996) are listed in Table 2. The moment of EGF3, which was not determined by the F-net, was regarded as the same as that of EGF1, because the averaged spectral ratio of EGF3 with respect to EGF1 at low frequencies at the target stations was almost unity (left panel of Fig. 3). Similarly, the moment of EGF4, which was again not determined by the F-net, was regarded as the same as that of EGF2, because the averaged spectral ratio of EGF4 with respect to EGF2 at low frequencies at the target stations was almost unity (right panel of Fig. 3).

Fig. 3

Fourier spectral ratios between the EGF events. The left panel shows the spectral ratio of EGF3 with respect to EGF1 at eight target stations (gray lines) and their logarithmic mean (green line). The mean spectral ratio is almost unity at low frequencies. Similarly, the right panel shows the spectral ratio of EGF4 with respect to EGF2 at eight target stations (gray lines) and their logarithmic mean (green line). The mean spectral ratio is almost unity at low frequencies

The corner frequencies of the EGF events were determined based on the observed Fourier spectral ratios between the main shock and the EGF events (Fig. 4). The gray lines in each panel show the spectral ratios of the main shock with respect to the EGF event at eight target stations and the green line displays their logarithmic mean. The following function based on the omega-square model was fitted to the mean spectral ratio:

where M01 is the moment of the large event, M02 is the moment of the small event, fc1 is the corner frequency of the large event, and fc2 is the corner frequency of the small event. The above-mentioned values of the seismic moment were used. The red line in each panel represents the fitted curve. The corner frequencies were evaluated as 3.9, 2.8, 2.7, and 3.0 Hz for EGF1–EGF4, respectively. Thus, we ascertained that the target frequencies of this study are lower than the corner frequencies of the EGF events.

Fig. 4

Fourier spectral ratios between the main shock and EGF events. The gray lines in each panel display the spectral ratios of the main shock with respect to the EGF event at eight target stations, the green line shows their logarithmic mean, and the red line is the fitted curve assuming the omega-square model

The inversion analysis in Case 1 was repeated for different values of the first-time-window triggering velocity. The goodness of fit was then evaluated by introducing variance reduction (VR), defined by:

where s(t) and o(t) denote the synthetic and observed waves, respectively. The integration was performed on the same time sections as used for the inversion (Fig. 2). Figure 5 shows the relation between the first-time-window triggering velocity and VR; VR takes the maximum value at the triggering velocity of 2.1 km/s. The results corresponding to this triggering velocity will be presented in the following sections.

Fig. 5

Relation between the first-time-window triggering velocity and the variance reduction VR (%). The VR is defined in Eq. (8). The results correspond to Case 1; VR takes the maximum value at the triggering velocity of 2.1 km/s

The synthetic velocities (0.2–2 Hz) were compared with those observed in Case 1 (Fig. 2). The agreement between the observed and synthetic waves is satisfactory at many stations. At KMMH16, the coincidence of the NS components is almost perfect including the pulses, while the amplitude of the EW component is underestimated, which is probably due to the fact that the coincidence of the mechanism between the large and the small evens is not perfect. Similarly, the coincidence of the NS component is almost perfect at KMM011, while the EW component is underestimated. The result for the EW component at KMM006 is good, while the NS component is underestimated.

The synthetic velocities (0.2–2 Hz) for Cases 2, 3, and 4 were compared with the observed ones (Figs. 6, 7, 8). The agreement between the observed and synthetic waves is satisfactory. In Cases 2 and 4, the result for the EW component at KMM011 was improved when EGF3 was used for the western part of the fault, while the result for the NS component at the same station was degraded. The results indicate that the fault mechanisms of both EGF1 and EGF3 do not perfectly coincide with the western part of the main shock. In Cases 3 and 4, the result for KMM004 was degraded when EGF4 was used for the eastern part of the fault, while the result for KMMH06 was improved. The variance reduction was 62.4, 65.8, 61.4, and 63.9% for the Cases 1, 2, 3, and 4, respectively. The largest variance reduction was obtained when EGF3 and EGF2 were combined.

Fig. 6

Comparison of the observed and synthetic velocity waveforms (0.2–2 Hz) for Case 2. The gray and black traces are the observed waveforms. The black traces correspond to the 15-s time sections used in the waveform inversion. The red traces are the synthetic waveforms

Fig. 7

Comparison of the observed and synthetic velocity waveforms (0.2–2 Hz) for Case 3. The gray and black traces are the observed waveforms. The black traces correspond to the 15-s time sections used in the waveform inversion. The red traces are the synthetic waveforms

Fig. 8

Comparison of the observed and synthetic velocity waveforms (0.2–2 Hz) for Case 4. The gray and black traces are the observed waveforms. The black traces correspond to the 15-s time sections used in the waveform inversion. The red traces are the synthetic waveforms

Figure 9 shows the final slip distribution for the 2016 Kumamoto earthquake obtained from the waveform inversion Cases 1–4. The star indicates the JMA hypocenter, which is the rupture starting point assumed for the waveform inversion. The absolute value of the slip was derived by using the moments of the EGF events and assuming a shear wave velocity of 3.55 km/s and density of 2.4 g/cm3. These values correspond to the second layer of the subsurface structure employed by Fukuyama et al. (1998). All slip models correspond to Mw7.2. Note that the dependence of the final slip distribution on the selection of the EGF events is small. The results clearly indicate that a significantly large slip region was located approximately 15 km northeast of the hypocenter. This location is similar for the different cases.

Fig. 9

Final slip distribution for the 2016 Kumamoto earthquake for Cases 1–4. The star indicates the JMA hypocenter, which is the rupture starting point assumed for the waveform inversion. The absolute value of the slip was obtained by using the moments of the EGF events. All of the results correspond to Mw7.2

Figure 10 shows the peak slip velocity distribution for the 2016 Kumamoto earthquake for Cases 1–4. The slip velocity for each subfault for each time window was estimated by dividing the amount of slip with the width of the time window (0.25 s). The largest velocity was then selected from the 12 time windows. Thus, it should be noted that the peak slip velocity shown here is a value averaged over the time window of 0.25 s. Again, note that the dependence of the peak slip velocity distribution on the selection of the EGF events is small. Figure 10 clearly shows that the slip velocity was largest approximately 15 km northeast of the hypocenter, where the final slip was also largest.

Fig. 10

Peak slip velocity distribution for the 2016 Kumamoto earthquake for Cases 1–4. The star indicates the JMA hypocenter, which is the rupture starting point assumed for the waveform inversion. The slip velocity of each subfault for each time window was estimated by dividing the amount of slip by the width of the time window (0.25 s). The largest velocity was then selected from the 12 time windows. Thus, it should be noted that the peak slip velocity shown here is a value averaged over the time window of 0.25 s

The inversion results reported in the previous section clearly indicate that a region of significantly large slip and slip velocity existed approximately 15 km northeast of the hypocenter. This will be referred to as “asperity” in the following discussion. All four slip models obtained from different combinations of EGFs suggest that the rupture, which was initiated at the hypocenter, first proceeded to the deeper part (red arrow in Fig. 11). The rupture then changed its direction and proceeded to the shallower part, across the “asperity” (Fig. 11). This process can be confirmed in Fig. 12, which shows the spatiotemporal distribution of the slip for Case 1. Although the slip and slip velocity in the western part of the fault plane depend more or less on the selection of the EGFs, in the region between the hypocenter and KMMH16, both the slip and slip velocity were relatively small; no “asperity” was observed between the hypocenter and KMMH16. Thus, it is not appropriate to conclude that the large-amplitude pulse-like ground motion in Mashiki was generated by the forward-directivity effect associated with the rupture of an asperity; the generation mechanism was different from that of the Kobe ground motion. As far as the source effect is concerned, the ground motion in Mashiki cannot be interpreted as the worst case scenario. On the other hand, the rupture of the “asperity” 15 km northeast of the hypocenter should have caused significantly large ground motions in regions close to the asperity. The significant damage of highway bridges in the region, including the Ogiribata Bridge (Nikkei Construction 2016; Fig. 1), can potentially be attributed to the inertia force caused by the rupture of the asperity, although strong ground motions also depend on the site effects.

Fig. 11

Peak slip velocity distribution for the 2016 Kumamoto earthquake for Case 1. The red arrows schematically illustrate the progression of the rupture

Fig. 12

Spatiotemporal distribution of slip for Case 1. The slip is shown for every 2 s

2.

Potential effects of soil nonlinearity on the waveform inversion results

As mentioned in “Background,” near-source stations are used in the present analysis. The inversion results could therefore be affected by soil nonlinearity because the amplitudes of strong ground motions were large in the near-source region of the Kumamoto earthquake. As a simple test of the effect of soil nonlinearity on the inversion results, velocity waveforms at relatively distant stations, which were not used in the inversion analysis, were simulated using the inversion results. The idea is that if the slip distributions are not contaminated by soil nonlinearity, they could be applied to distant stations, where the ground motions are less affected by soil nonlinearity. The same combination of EGF events utilized in the Case 1 inversion and the stations shown in Fig. 1 (open triangles) were used for the Case 1 slip model. The simulations were performed in the same frequency range (0.2–2 Hz). Borehole records were used for the KiK-net stations. The results were compared with the observed waveforms (Fig. 13). The simulation results explain the main features of the observed waveforms well, except for the NS component at KMMH03, if we consider the fact that the stations shown in Fig. 1 were not used in the inversion. However, one could point out the arrival time delay of the synthetic waveforms with respect to the observed ones. The delay is especially significant for the latter part of the waveforms. These interesting phenomena could be interpreted as follows. In general, soil nonlinearity accompanies the reduction in the shear wave velocity within the soil. This will lead to the delay of the arrival times of seismic waves (Nozu and Morikawa 2003). The rupture time could therefore be delayed in the inversion results accounting for the late arrival times if near-source strong-motion records are contaminated by soil nonlinearity. Nevertheless, Fig. 13 shows that the delay of the rupture time does not exceed 1 s. The variance reduction in the inversion analysis was maximized when the first-time-window triggering velocity of 2.1 km/s was used. One could argue that the velocity was slightly smaller than expected for a crustal earthquake. It might be possible to assign a larger value if we could eliminate the effects of soil nonlinearity. In addition, because of this potential soil nonlinearity effect, the authors do not insist on the moment magnitude obtained in this study (Mw7.2), which is slightly greater than the results of other studies such as the F-net moment tensor solution and Asano and Iwata (2016). The soil nonlinearity effect could lead to an increase in the low-frequency level of the observed ground motions, which could cause the overestimation of the seismic moment.

Fig. 13

Comparison of the observed and synthetic velocity waveforms (0.2–2 Hz) at the distant stations. The locations of the stations are shown in Fig. 1. These stations were not used in the inversion analysis. The black traces are the observed waveforms. The red traces are the synthetic waveforms. The Case 1 slip model was used

3.

Comparison with other source studies

As mentioned in “Background,” the use of empirical Green’s functions for the analysis of the Kumamoto earthquake is quite advantageous if we consider the 3D nature of the subsurface structure in the near-source region of the Kumamoto earthquake represented by the existence of the Aso caldera. On the other hand, there are potential disadvantages in using empirical Green’s functions for the inversion analysis. Table 3 summarizes the potential sources of errors for each of the numerical and empirical Green’s functions. The primary source of errors for the numerical Green’s functions is the uncertainty in the subsurface structure model. To mitigate this problem, subsurface structure model tuning using records of small events is often employed (e.g., Hikima and Koketsu 2005; Asano and Iwata 2011). The primary source of errors for the empirical Green’s functions is associated with the allocation of small events to the subfaults. As mentioned in “Methods” section, to assign empirical Green’s functions to each of the subfaults, conventional corrections for the geometrical spreading and time shifts (Irikura 1983) are applied to the empirical Green’s functions in Eq. (2). However, there are some aspects that cannot be accommodated by this type of correction. For example, it is not easy to incorporate the frequency-dependent radiation coefficient in the allocation of empirical Green’s functions. This problem can be mitigated to some extent by dividing the fault plane into two or more parts and by assigning different small events to those parts (e.g., Nozu 2007; Nozu and Irikura 2008). In summary, each Green’s function has its own source of errors and mitigation method. However, it will never be possible to obtain “perfect” Green’s functions. With respect to this fact, the diversity of the inversion method will be one of the keys to ensure the reliability of the inversion results. Therefore, it is important to compare the inversion results obtained in this study using the empirical Green’s functions to those obtained utilizing different types of Green’s functions. Hence, the result of this study was compared with an inversion result obtained from numerical Green’s functions for a layered half-space that was appropriately tuned for each site (Asano and Iwata 2016). The two results share the main features in spite of the different Green’s functions and stations used. For example, a significantly large slip region was located approximately 15–20 km northeast of the hypocenter in the model of Asano and Iwata (2016), which is fairly close to the large slip region of the present study (Fig. 9). Therefore, it can be concluded that these two source models capture the main features of the rupture process of the earthquake.

Table 3

Potential sources of errors for each of the numerical and empirical Green’s functions

Green’s functions

Potential sources of errors

Mitigation method

Numerical Green’s functions

Uncertainty in the subsurface structure model

Tuning of the subsurface structure model using records of small events (e.g., Hikima and Koketsu 2005; Asano and Iwata 2011)

Empirical Green’s functions

Errors associated with the allocation of small events to the subfaults

The fault plane is divided into two or more parts and different small events are assigned to those parts (e.g., Nozu 2007; Nozu and Irikura 2008)

On the other hand, one important difference between the two models is that the model of Asano and Iwata (2016) indicates a shallow slip near the station KMMH16. In contrast, such a shallow slip patch was not needed in our study to fit the waveform at KMMH16. Asano and Iwata (2016) included low-frequency components (0.05–0.2 Hz) in their analysis, which were not used in the present analysis. The velocity waveforms in this frequency range observed at KMMH16 are closely related to a fling step. In fact, the velocity waveforms, which were successfully reproduced by the model of Asano and Iwata (2016), clearly exhibit a one-sided nature. Therefore, the shallow slip patch appearing in their model was presumably needed to explain the fling step. On the other hand, the present source model includes information related to the 1-Hz energy observed at KMMH16.

4.

Effects of the dip angle

One of the main differences between the present source model and those of other studies could be the dip angle. The present model uses the dip angle of 84°, while some of the other source models, including that of Asano and Iwata (2016), employ a lower dip angle. Therefore, the effect of the dip angle was investigated.

A fault plane with a low dip angle that involves the JMA hypocenter does not match the surface fault trace as long as the strike angle is fixed to 232°. Therefore, using a low dip angle inevitably requires the use of two fault planes with different strike angles. Then, as shown in the top two panels of Fig. 14, in addition to the fault plane with a strike angle of 232°, a fault plane with strike and dip angles of 205° and 72°, respectively, was assumed (Asano and Iwata 2016). This fault was extended from the hypocenter toward the northeast, for 2 km (top left) or 4 km (top right), and then connected to the fault with the strike angle of 232°. In the top-left case, the latter fault had a dip angle of 78°, consistent with the surface fault trace. In the top-right case, the latter fault had a dip angle of 72°. The other inversion parameters, including the number of subfaults, are exactly the same as in Case 1. A triggering velocity of 2.1 km/s was used.

Fig. 14

Additional inversion cases to examine the effect of the strike and dip angles. The top two panels show the fault models with two fault planes, one having a strike angle of 205° and the other having a strike angle of 232°. The former fault was extended from the hypocenter toward the northeast, for 2 km (top left) or 4 km (top right), and then connected to the fault with the strike angle of 232°. In the top-left case, the latter fault had a dip angle of 78°, consistent with the surface fault trace. In the top-right case, the latter fault had a dip angle of 72°. The bottom two panels show the corresponding final slip distributions

The bottom two panels of Fig. 14 show the corresponding final slip distributions. The results are similar to that of Case 1 in spite of the different strike and dip angles; the results imply that the basic features of the slip models (such as the large slip approximately 15 km northeast of the hypocenter) do not significantly depend on the strike and dip angles. The variance reduction was 60.6 and 60.3% for the dip angle of 78° and 72°, respectively. The results are slightly inferior to that of Case 1 with a variance reduction of 62.4%. However, because the difference is small, the authors do not insist that a nearly vertical fault plane is more appropriate.

The rupture process of the main shock of the 2016 Kumamoto earthquake was investigated in this study based on the inversion of strong ground motions, particularly the generation of strong ground motions in the frequency range relevant to structural damage. Strong-motion records in the near-source region were mainly used because the authors were interested in the generation mechanism of the damaging ground motions in the near-source region. Empirical Green’s functions (EGFs) were employed to avoid uncertainty in the subsurface structure model. Four cases of inversions with different combinations of small events were analyzed to investigate the dependence of the inversion results on the selection of the small events. It was found that the dependence of the final slip and peak slip velocity distributions on the selection of the EGF events was small. The results clearly indicate that a region of significantly large slip and slip velocity existed approximately 15 km northeast of the hypocenter. No “asperity” was observed between the hypocenter and Mashiki. Thus, it is not appropriate to conclude that the large-amplitude pulse-like ground motion in Mashiki was generated by the forward-directivity effect associated with the rupture of an asperity. As far as the source effect is concerned, the ground motion in Mashiki cannot be interpreted as the worst case scenario. On the other hand, the rupture of the “asperity” 15 km northeast of the hypocenter should have caused significantly large ground motions in regions close to the asperity. The significant damage of highway bridges in the region, including the Ogiribata Bridge, can potentially be attributed to the rupture of the asperity. The results of this study were compared with an inversion result obtained using numerical Green’s functions for a layered half-space. The two results share the main features in spite of the different Green’s functions and stations used. Therefore, it can be concluded that these two source models capture the main features of the rupture process of the earthquake.

Authors’ contributions

AN conducted the waveform inversion. YN participated in the discussion on the relation between the rupture process and the near-source ground motions. Both authors read and approved the final manuscript.

Acknowledgements

The authors are grateful to the National Research Institute for Earthquake Science and Disaster Resilience, Japan, for providing strong-motion data of K-NET and KiK-net and the moment tensor solutions by the F-net. The authors appreciate valuable comments from Dr. Haruo Horikawa, Dr. Arben Pitarka, and an anonymous reviewer.

Competing interests

The authors declare that they have no competing interests.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (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.