# Bulletin of the Seismological Society of America

## Abstract

The *M*_{w} 7.0 Kumamoto, Japan, earthquake occurred on 15 April 2016 at 16:25 UTC. Using ground accelerations recorded by 104 near‐field stations, we investigate spatial variability of observed ground motions, apparent period dependence, and azimuthal variation, as well as rupture directivity effects on various intensity measures. We develop a simplified ground‐motion model that includes both geometric and anelastic attenuation terms. Comparisons of observed and predicted ground motions suggest that predictions from the Next Generation Attenuation‐West2 models provide good fits for the overall observation. Analysis of spatial distribution of the residuals shows that observed peak ground velocity (PGV) and long‐period spectral accelerations (SAs) in the 150°–180° azimuth range along the rupture backward direction (southwest of the fault) can be as low as 0.3–0.8 times the average observation of this event. Long‐period ground motions on the northeast side of the fault in the forward direction are much higher than average, with PGV and long‐period SAs ranging from 1.2 to 1.5 times the average. There is clear period dependence of the strong ground motion variation. The biases due to directivity generally decrease with decreasing period for all azimuth ranges. On the distance dependence of directivity effects, our study shows that directivity effects can be considered practically nonsignificant for stations close to the hypocenter. We also perform a log–linear regression of the residuals, using a new directivity predictor. Our results show that for the 2016 *M*_{w} 7.0 Kumamoto earthquake, rupture directivity produces significant amplifications in the rupture forward direction, whereas deamplification effects are observed in the rupture backward region. Directivity effects are particularly relevant for PGV and long‐period SA (i.e., SA at periods ≥2.0 s). Such effects do not have systematic influence on peak ground acceleration and short‐period ground motions (i.e., SA at periods <2.0 s).

*Electronic Supplement:*Figures of variation of regression residuals with *R*_{rup} for observed peak ground acceleration (PGA), peak ground velocity (PGV), and spectral accelerations (SAs) and of regression residuals versus *V*_{S30}.

## Introduction

On 15 April 2016, 16:25 UTC, an earthquake of magnitude *M*_{w} 7.0 struck Higashi Ward in the city of Kumamoto (2010 population: ∼731,000) in the Kyushu region in southwest Japan. Significant damages were experienced in the following prefectures: Kumamoto, Oita, Fukuoka, Yamaguchi, Saga, Nagasaki, and Miyazaki. The prefecture of Kumamoto experienced the worst consequences, with over 50,000 buildings damaged. Among them, more than 2400 totally collapsed, according to the prefectural government. Various seismically induced ground failure mechanisms occurred in this prefecture (i.e., landslides, liquefaction, surface fault rupture); also, fire following the earthquake was reported. Important lifelines were affected by the earthquake, including power, water, and gas systems. At least 66 fatalities and over 1500 injuries were reported. Early estimates of the economic costs of the damage ranged from $5.5 billion U.S. to $7.5 billion U.S., with insured property losses estimated to be between $800 million U.S. and $1.2 billion U.S.

The *M*_{w} 7.0 Kumamoto earthquake occurred at a depth of ∼10.0 km in the Futagawa and Hinagu fault zones. Twenty‐eight hours before the mainshock, on 14 April, an *M*_{w} 6.1 foreshock had shaken the same region severely. Koji (2016) shows that the Futagawa fault was the source of the *M*_{w} 7.0 earthquake, with 2 m+ right‐lateral strike slip at the surface; the *M*_{w} 6.1 foreshock ruptured for about 15 km in the northernmost section of the Hinagu fault, without surface rupture evidences. Inverse analysis of the rupture process suggests that the *M*_{w} 7.0 Kumamoto earthquake had a strong northeast rupture directivity (Asano and Iwata, 2016; U.S. Geological Survey [USGS], 2016; Yagi *et al.*, 2016). These findings are consistent with the observed distribution of strong ground shaking and with the orientation of the mapped faults in the region.

The Kumamoto mainshock was well recorded by K‐NET and KiK‐net, two strong‐motion networks in Japan operated by the National Research Institute for Earth Science and Disaster Resilience (NIED). According to NIED, over 690 groups of three‐component (east–west, north–south, and up–down) accelerograms were recorded by triggered strong‐motion stations. A large number of digital strong‐motion recordings were available in the near‐fault region. Over 75 groups of accelerograms were recorded by stations located within 60 km of the Futagawa fault. The largest peak ground acceleration (PGA) was 1157 cm/s^{2} in the east–west direction. It was recorded by the KMMH16 station (Fig. 1), close to the epicenter. These near‐fault recordings provide an invaluable opportunity for studying rupture directivity effects.

During the Next Generation Attenuation‐West2 (NGA‐West2) project (Bozorgnia *et al.*, 2014), five directivity models were developed. These models are based on both data from the NGA‐West2 database and numerical simulations. However, Spudich *et al.* (2013, 2014) show that a large discrepancy still exists among these models, especially for dipping faults. As a result, directivity models are not incorporated into the NGA‐West2 ground‐motion models (GMMs), with the exception of the Chiou and Youngs (2014) model.

In this article, a total of 104 groups of three‐component acceleration recordings are used to investigate the rupture directivity effects on strong motion during the *M*_{w} 7.0 Kumamoto earthquake. These records are examined and compared with estimations from the NGA‐West2 models, based on a global database (Ancheta *et al.*, 2014). We investigate spatial variability of the strong ground motion. We also analyze rupture directivity effects, such as systematic azimuthal variation effects and apparent period dependence of the near‐fault strong motion, using residual analysis. Based on these results, we present a quantitative directivity model for this event which includes a new geometric directivity predictor. We anticipate that results from this study will be useful for future GMMs development, as well as for a broader improvement of the technical knowledge on directivity effects for large shallow crustal earthquake events.

## Strong‐Motion Data

The strong‐motion dataset used in this study consists of three‐component accelerograms recorded by the NIED strong‐motion seismograph networks of Japan. These data are available in a unified website for both K‐NET and KiK‐net networks (see Data and Resources). Figure 1 shows the 104 stations used in this study. The selected stations are all located within 80 km of the finite‐fault plane of the *M*_{w} 7.0 Kumamoto earthquake. According to Japan Meteorological Agency (JMA), an *M*_{JMA} 5.7 event occurred ∼30 s after the mainshock (Yoshida, 2016) in the central part of the Oita prefecture (see Fig. 1). Observed ground motions in the Oita prefecture show waveforms with dominant high‐frequency components after the mainshock *S*‐wave arrival. These high‐frequency motions are related to the *M*_{JMA} 5.7 aftershock (Aitaro *et al.*, 2016; Yoshida, 2016). As a result, we checked all time series affected by this aftershock. We identified the mainshock and the *M*_{JMA} 5.7 aftershock using the *P*‐ and *S*‐wave arrival in the acceleration and velocity time series for each recording in the Oita prefecture (station code OIT). Figure 2a,b shows the arrival time of the *M*_{JMA} 5.7 aftershock for the acceleration and velocity of the east–west component waveforms, respectively.

For KiK‐net stations, we only use three‐component recordings observed on the ground surface, because intensity measures (IMs) recorded by downhole accelerometers are usually lower than those recorded on the surface, and GMMs are developed using data observed on ground surface from free‐field stations. All selected recordings from K‐NET and KiK‐net are sampled at 100 samples per second. The K‐NET and KiK‐net unified website also provides soil profile data at strong‐motion station locations. The sites are classified based on the time‐averaged shear‐wave velocity in the upper 30 m (*V*_{S30}). For K‐NET stations with soil profile depth less than 30 m, we use a *V*_{SZ}‐to‐*V*_{S30} extrapolation model for Japan (Boore *et al.*, 2011) to calculate *V*_{S30}, in which *V*_{SZ} is the average shear‐wave velocity of soil layers from surface to depth *Z*. Source‐to‐site distances for each station are calculated using the earthquake source geometry (Fig. 1) provided by Yagi *et al.* (2016). To avoid missing displacement information after high‐pass filtering and recover reliable permanent displacements, we process the records using a modified version of the multi‐consecutive‐segment baseline correction method (e.g., Boore, 2001; Akkar and Boore, 2009), initially proposed by Iwan *et al.* (1985).

## Comparison of Observed Horizontal Ground Motion with NGA‐West2 Ground‐Motion Models

The IMs used in this study are PGA, peak ground velocity (PGV), and spectral acceleration (SA) at various periods. We define them as the geometric mean of PGAs, PGVs, and SAs of the two as‐recorded horizontal (east–west and north–south) components of the ground motion. In this section, we compare observed variation of the horizontal ground‐motion IMs with distance, with predictions from the NGA‐West2 GMMs. We perform these analyses using the closest distance to the rupture plane *R*_{rup} as source‐to‐site distance. The observed data are adjusted to a common velocity using the Seyhan and Stewart (2014) model (hereafter, SS14, which represents the site term in the Boore *et al.*, 2014, model), in which amplitude and frequency‐dependent site effects are a function of *V*_{S30}, relative to a reference site condition (i.e., 760 m/s).

Based on strong‐motion data introduced in the Strong‐Motion Data section, we develop an *ad hoc* simplified GMM that includes both geometric and anelastic attenuation terms, with rupture distance *R*_{rup}. The functional form, similar to the Campbell and Bozorgnia (2014) model, is given as (1)in which *Y* represents the observed IM, and coefficients *a*, *b*, *c*, and *d* are derived from regression analysis performed using the least‐squares method; *σ*_{ln(Y)} is the standard deviation of the regression.

Table 1 shows regression analysis results for PGA, PGV, and SA for the following periods: 0.2, 0.5, 1.0, 2.0, 3.0, 5.0, 7.5, and 10.0 s. Regression coefficients, determination coefficients (*R*^{2}), and standard deviation of residuals are shown in Table 1 as a function of the considered IM. The determination coefficients *R*^{2} provide a measure of how well observed data are fitted by the regression model. Ⓔ Figure S1 (available in the electronic supplement to this article) shows regression residuals versus *R*_{rup}, indicating that our regression results fit the data well. Ⓔ Figures S2 and S3 show regression residuals versus *V*_{S30} for PGA, PGV, and SA for the following periods: 0.2, 1.0, 3.0, and 5.0 s, before and after the site term adjustment (i.e., data adjusted to a common velocity using the SS14 model). Ⓔ Figure S2 shows a trend of decreasing residuals with increasing velocity. The impact of the site term adjustment (SS14 model) is clear in Ⓔ Figure S3, which shows residuals roughly aligned around zero. Figure 3 shows variation of the considered IMs (PGA, PGV, and SA) with *R*_{rup}, along with least‐squares regression lines obtained from equation (1). These results are compared with median predictions from the NGA‐West2 models (Abrahamson *et al.*, 2014, hereafter, ASK14; Boore *et al.*, 2014, hereafter, BSSA14; Campbell and Bozorgnia, 2014, hereafter, CB14; Chiou and Youngs, 2014, hereafter, CY14). The predictions from the NGA‐West2 models provide a reasonably good fit for the observed data, with the exception of PGV, for which predictions are consistently lower than observations.

Figures 4 and 5 show observed acceleration spectra in the forward and backward direction of fault rupture respectively, along with predicted SAs versus period (*T*) from NGA‐West2 models. Location of strong‐motion stations along with their codes are shown in Figure 1. Figure 4 shows that, in the forward direction, mean predictions from NGA‐West2 models underestimate observed spectra for long‐period ground motions (i.e., *T*≥2.0 s) at the following stations: KMM004, KMM005, KMMH06, OIT006, OIT008, OIT010, and OIT012. Observed spectrum at station KMMH02 is uniformly higher than mean predictions from NGA‐West2 models for all periods (0.01–10 s). These results suggest that, in the forward direction, rupture directivity increases ground shaking intensity levels, especially, but not only, for long periods. For stations KMM006 and KMMH16, located close to the hypocenter in the forward direction, NGA‐West2 models provide relatively good predictions at long periods, indicating that rupture directivity effects are nonsignificant at these two locations. For recordings affected by the triggered *M*_{JMA} 5.7 aftershock (e.g., OIT006, OIT010, and OITH05), we compared the calculated spectra before and after removing the aftershock waveforms. Our results show that the aftershock only affects spectra at short periods (i.e., *T*<2.0 s). For long periods (i.e., *T*≥2.0 s), the difference is practically negligible. These observations confirm that amplification effects on long‐period ground motions are not affected by the aftershock.

Figure 5 shows that in the backward direction, NGA‐West2 model predictions overestimate observed spectra at stations KMM019 and KMM021 for almost all periods. Observed spectra at the following stations: KMM018, KMM020, KMM022, KMMH10, KGS001, and NGS014 are generally lower than those predicted by the NGA‐West2 models, for long‐period ground motions (*T*≥2.0 s).

## Spatial Variability of Observed Ground Motions

Previous studies have shown that near‐fault ground motion is significantly influenced by the focal mechanism and rupture process of causative fault (e.g., Somerville *et al.*, 1997; Bradley and Cubrinovski, 2011; Xie *et al.*, 2014). In this study, we use residual analysis for investigating the spatial distribution, the apparent period dependence of observed near‐fault strong motion, and to relate directivity effects to azimuthal variations. The residuals are calculated as the difference between each observed ground‐motion IM, and the estimation performed using the mean regression relationship for this event (equation 1) is (2)in which *Y*_{obs} is the observed IM adjusted to *V*_{S30}=760 m/s, using the SS14 model for each site, and is the mean predicted value of each IM determined using equation (1).

Figure 6 shows the spatial distribution of residuals for various IMs. The smallest residuals (−2.0 to −1.5) are located southwest of the causative fault (rupture backward direction). Thus, ground motions observed on the opposite side of the rupture direction are significantly lower than the average values. Ground motions on the northeast side of the fault (rupture forward direction) are generally higher than average, especially for long periods, with largest residuals equal to 2.0. In particular, the residuals of PGV and SA at 3.0 and 5.0 s exhibit clear directivity amplification in the rupture forward region, whereas residuals of PGA and short periods do not. The observed characteristics of residuals spatial distribution are consistent with those observed during the 2008 *M*_{w} 7.9 Wenchuan event and the 2014 South Napa earthquake (Wen *et al.*, 2010; Baltay and Boatwright, 2015).

Another way to investigate the spatial variability of ground motions is to plot acceleration and velocity waveforms against distance for near‐fault stations (Fig. 2). In the forward directivity region, both acceleration and velocity increase very rapidly. Acceleration and velocity waveforms are characterized by higher peaks and shorter durations in the forward region. Velocity pulses are clear in the forward direction. Results from this analysis show evident effects of rupture directivity and its dependence on rupture distance.

## Quantification of Rupture Directivity Effects, Apparent Period Dependence, and Azimuthal Variation Effects

Spudich *et al.* (2013) show that the use of normalized site distances produces undesirable effects on geometry definition, especially for relatively large earthquakes associated with large rupture dimensions. As a result, differently from previous models (e.g., Somerville *et al.*, 1997), the directivity models developed during the NGA-West2 project use source‐to‐site distances (e.g., rupture distance *R*_{rup} or Joyner–Boore distance *R*_{JB}) in kilometers (e.g., distances are not normalized to fault length). In this study, we use the same geometry definitions used in the Bayless and Somerville (2013) model. However, based on observation during this event, we define a different geometric directivity predictor (*f*_{g}): (3)in which *s* is the along‐strike projected length of the rupture surface toward the site (in km); and *θ* is the angle (°) between rupture direction and the direction of source to sites. In the definition of *θ*, we use direction of dominant rupture propagation (northeast) as rupture direction for this strike‐slip event. Equation (3) is defined for *θ* between 0° and 180°. Variable *s* is set to be 1 km for sites with rupture distance *s* <1 km, so that the term ln*s* has a minimum value of zero. Our geometric directivity predictor (equation 3) is defined using observation within 80 km from the finite fault of the *M*_{w} 7.0 Kumamoto earthquake. It captures both amplification effects in the forward directivity direction and deamplification effects in the backward directivity direction. In both cases, using our geometric predictor, the strength of directivity effects increases as *s* increases, with a logarithmic trend (i.e., for a fixed value of *θ*, its rate decreases with increasing *s*). We recommend applying a distance taper to equation (3). Bayless and Somerville (2013) provide some guidance on the choice of this value in the forward directivity direction (discussed later in this section). We recognize that the use of a broad, global dataset with a large number of events is needed to further verify these observations, and to define a proper distance taper, also in the backward directivity zone.

Directivity effects are inherently related to systematic azimuthal variations in the strong‐motion amplitude residuals. Figure 7 shows residuals versus azimuth *θ* of the observation station sites for various IMs. For PGV and long‐period SAs at 2.0, 3.0, 5.0, 7.5, and 10.0 s, there is a clear trend of decreasing residuals with increasing azimuths in the 0°–180° range (Fig. 7b). This systematical trend is not present for PGA and short‐period SAs at 0.2, 0.5, and 1.0 s (Fig. 7a).

Based on the values of azimuth *θ*, we group residual data into six bins. Table 2 shows mean residuals and standard deviations for these azimuth bins. In the 120°≤*θ*<180° range, the mean residuals are uniformly negative for both long‐period and short‐period motion. In the 0°≤*θ*<60° range, residuals are uniformly positive for PGV and long‐period SA at 3.0, 5.0, 7.5, and 10.0 s. Residuals for short‐period motions in this azimuth range do not have a clear and consistent trend. Observed ground motions of PGV and long‐period SA at 3.0, 5.0, 7.5, and 10.0 s are as low as 0.3–0.8 times the mean fit of this event in the 150°≤*θ*<180° range, and equal to 1.2–1.5 times the mean in the 0°≤*θ*<30° range. These results suggest that rupture directivity has a significant effect on the spatial distribution of near‐fault strong motions for this event.

We believe that the evaluation of a possible period dependence of directivity effects is critical for ground-motion characterization and seismic-hazard assessments. As a result, we perform an accurate analysis of this feature based on visual inspection of residuals and side‐to‐side comparisons of various azimuth groups. Figure 8 shows variation of mean residuals with period for various azimuth groups. There is clear period dependence of the strong ground motion variation. In particular, there is a trend of increasing residuals with increasing period in the forward direction (e.g., azimuth range 0°≤*θ*<60°). Conversely, the biases due to directivity generally decrease, or remain flat, with decreasing period for all azimuth ranges. These trends show that, for this event, the strong ground motion has an amplification in the rupture forward direction for long periods, whereas there is weakening of the directivity effect at short periods. These observations are consistent with recent findings by Gallovič (2016) for the 2014 *M*_{w} 6.0 South Napa earthquake and Pacor *et al.* (2016) for weak‐motion data from Italy.

Ⓔ Figures S4 and S5 show residuals versus rupture distance parameter *s* for PGV and long‐period SA at 2.0, 3.0, 5.0, 7.5, and 10.0 s, for *R*_{rup}<*L*, and *R*_{rup}>*L*, respectively (in which *L* is equal to 56 km, the length of the ruptured fault). Ⓔ Figures S4 and S5 provide a confirmation that directivity effects are nonsignificant (residuals are close to zero) for stations located close to the hypocenter (i.e., |*Ln*(*s*)|<1 in Ⓔ Figs. S4 and S5), although this effect is not clear for all ground‐motion IMs. Furthermore, residuals in the forward direction increase with rupture distance, whereas residuals in the backward direction decrease with increasing distance, for PGV and long‐period SA at 2.0, 3.0, 5.0, 7.5, and 10.0 s. Ⓔ Figure S5 also shows that directivity effects can still be strong when *R*_{rup}>*L*, especially for PGV and long‐period SA at 7.5 and 10.0 s. This suggests that the distance taper function of the Bayless and Somerville (2013) model may not be appropriate for this event.

Figure 9 shows variation of residuals relative to cos(*θ*) for various IMs. It shows that residuals of PGV and long‐period SA at 2.0, 3.0, 5.0, 7.5, and 10.0 s are generally negative when *θ* is close to 180° (i.e., cos(*θ*)=−1, backward directivity zone) and positive when *θ* is close to 0 (i.e., cos(*θ*)=1, forward directivity zone). The residuals systematically increase with cos(*θ*) for PGV and long‐period SA at 2.0, 3.0, 5.0, 7.5, and 10.0 s; no such systematical trend is observed for PGA and short‐period SA at 0.2, 0.5, and 1.0 s. We fit the residuals using the following log–linear model which has been used by Bayless and Somerville (2013): (4)in which *f*_{D} is the residual of the natural logarithm of IM; *C*_{0} and *C*_{1} are the regression coefficients of the model. This model is valid for the geometric mean of the two as‐recorded horizontal (east–west and north–south) components of the ground motion. Figure 10 shows variation of residuals with the geometric directivity predictor (*f*_{g}) for PGA and short‐period SA at 0.2, 0.5, and 1.0 s (Fig. 10a), and PGV and long‐period SA at 2.0, 3.0, 5.0, 7.5, and 10.0 s (Fig. 10b). Our results show that only PGV and long‐period SA at 2.0, 3.0, 5.0, 7.5, and 10.0 s show a systematical trend of increasing residuals with increasing *f*_{g}. For PGV and long‐period SAs at 2.0, 3.0, 5.0, 7.5, and 10.0 s, mean regression curves of residuals versus the geometric directivity predictor (*f*_{g}) are also shown in Figure 10b. Regression coefficients, determination coefficients (*R*^{2}), and root mean square error (rmse) for PGV and long‐period SAs are shown in Table 3 as a function of the considered IM. We recommend using our results with particular care, because they are developed based on data from one specific event.

## Summary and Conclusions

The 2016 *M*_{w} 7.0 Kumamoto earthquake was well recorded by K‐NET and KiK‐net, strong‐motion networks in Japan. Over 104 groups of three‐component strong‐motion recordings are available in the near‐fault region of the earthquake. This substantial dataset provides invaluable data for the investigation of rupture directivity effects. Although several recent global directivity models are available in the literature, full consensus has not been reached among them (Spudich *et al.*, 2013, 2014).

In this study, we use the substantial dataset from the *M*_{w} 7.0 Kumamoto earthquake to investigate the following features: (1) spatial distribution characteristics of the ground motion and (2) rupture directivity effects on PGA, PGV, and SA at various periods. The observations are then compared with NGA‐West2 GMMs. We find that, for the whole dataset, NGA‐West2 models produce reasonably good estimation of the attenuation with distance for PGA, PGV, and SA at various periods. However, as a result of significant spatial variability induced by rupture directivity during this event, NGA‐West2 predictions underestimate long period (i.e., *T*≥2.0 s) SAs in the rupture forward direction by a large amount. In the backward direction, NGA‐West2 GMMs consistently overestimate observed spectra for all periods. For stations located close to the epicenter, NGA‐West2 models provide good predictions, indicating that rupture directivity effects can be considered nonsignificant.

Spatial distribution of the residuals shows that on the southwest side of the causative fault (rupture backward direction), observed ground motions are significantly lower than the average. On the northeast side of the fault (rupture forward direction), residuals are constantly higher than the average, with largest residual values equal to 2.0.

Residual analysis shows that observed PGV and long‐period SAs at 2.0, 3.0, 5.0, 7.5, and 10.0 s generally decrease with increasing values of azimuth *θ*. This systematical trend is not present for PGA and short‐period SAs at 0.2, 0.5, and 1.0 s. Residual values of long‐period SAs in the backward directivity zone with azimuth 150°≤*θ*<180° are as low as 0.3–0.8 times the mean fit of this event. In the forward directivity zone, long‐period ground motions can be as high as 1.2–1.5 times the mean in the 0°≤*θ*<30° range.

We observe an apparent period dependence of rupture directivity effects for this event. In particular, we show that, especially in the rupture forward direction, long‐period ground motions at periods equal to 2.0, 3.0, 5.0, 7.5, and 10.0 s are amplified, whereas for short periods (*T*<2.0 s), directivity effects are generally weakened. These findings are consistent with other recent studies carried out for earthquakes in California and Italy (i.e., Baltay and Boatwright, 2015; Gallovič, 2016; Pacor *et al.*, 2016). Our results show that effects of rupture directivity on long-period motions are azimuth‐ and distance-dependent. Directivity effects can be practically considered nonsignificant for locations close to the hypocenter. We then perform a log–linear regression of the residuals, developed using a new directivity predictor. The goodness of fit is then evaluated using confidence intervals, determination coefficients, and visual inspection of the data.

Even in recently developed directivity models, the rupture directivity effects are considered only for periods greater than 0.5 s. Thus, impact of rupture directivity on short periods is usually neglected (e.g., Bayless and Somerville, 2013). Our analyses provide a means by which to assess rupture directivity effect analyses for future GMMs and hazard analysis assessment.

## Data and Resources

The strong‐motion dataset consists of data recorded by the National Research Institute for Earth Science and Disaster Resilience (NIED) strong‐motion seismograph networks of Japan, available in a unified website for K‐NET and KiK‐net (http://www.kyoshin.bosai.go.jp/, last accessed June 2016). The soil condition data for strong‐motion stations are also obtained from the same website. The geophysical maps produced in this study are produced using Generic Mapping Tools (GMT; Wessel and Smith, 1991).

## Acknowledgments

We thank the National Research Institute for Earth Science and Disaster Resilience (NIED) in Japan for providing strong‐motion data for the *M*_{w} 7.0 Kumamoto earthquake. Financial support of this study was fully provided by the China National Special Fund for Earthquake Scientific Research (201408014) and Natural Science Foundation of China (51578513, 51208476, 51639006). Helpful comments from Associate Editor John Douglas, Frantisek Gallovič, and anonymous reviewers greatly improved the article. We also thank Jonathan P. Stewart for his help on the implementation of the Seyhan and Stewart (2014) model.

- Manuscript received 11 August 2016.