Kinematic Rupture Modeling of Ground Motion from the M7 Kumamoto, Japan Earthquake

We analyzed a kinematic earthquake rupture generator that combines the randomized spatial field approach of Graves and Pitarka (Bull Seismol Soc Am 106:2136–2153, 2016 ) (GP2016) with the multiple asperity characterization approach of Irikura and Miyake (Pure Appl Geophys 168:85–104, 2011 ) (IM2011, also known as Irikura recipe). The rupture generator uses a multi-scale hybrid approach that incorporates distinct features of both original approaches, such as small-scale stochastic rupture variability and depth-dependent scaling of rupture speed and slip rate, inherited from GP2016, and specification of discrete high slip rupture patches, inherited from IM2011. The performance of the proposed method is examined in simulations of broadband ground motion from the 2016 Kumamoto, Japan earthquake, as well as comparisons with ground motion prediction equations (GMPEs). We generated rupture models with multi-scale heterogeneity, including a hybrid one in which the slip is a combination of high- slip patches and stochastic small scale variations. We find that the ground motions simulated with these rupture models match the general characteristics of the recorded near-fault motion equally well, over a broad frequency range (0–10 Hz). Additionally, the simulated ground motion is in good agreement with the predictions from Ground Motion Prediction Equations (GMPEs). Nonetheless, due to sensitivity of the ground motion to the local fault rupture characteristics, the performance among the models at near-fault sites is slightly different, with the hybrid model producing a somewhat better fit to the recorded ground velocity waveforms. Sensitivity tests of simulated near-fault ground motion to variations in the prescribed kinematic rupture parameters show that average rupture speeds higher than the default value in GP2016 (average rupture speed = 80% of local shear wave speed), as well as slip rate durations shorter than the default value in GP2016 (rise time coefficient = 1.6), generate ground motions that are higher than the recorded ones at periods longer than 1 s. We found that these two parameters also affect the along strike and updip rupture directivity effects, as illustrated in comparisons with the Kumamoto observations.


Introduction
The characterization of earthquake rupture has been at the forefront of new developments in numerical methods for strong ground motion simulations of crustal and subduction zone earthquakes. Analyses of recent large earthquakes indicate that multi-scale earthquake rupture models that combine large-scale, large-slip features with smaller, but effective high-frequency seismic energy generation areas, better reproduce overall characteristics of ground motion on a broad frequency range Frankel et al. 2018;Iwaki et al. 2016a, b;Pulido et al. 2015;Dan et al. 2018;Miyake et al. 2003;Sekiguchi et al. 2008; Graves and Pitarka 2010;Pitarka et al. 2000). Schemes that include spatial correlations and interdependency of rupture parameters guided by rupture dynamics modeling offer significant advantages for use in broadband ground motion simulations (Imperatori and Mai 2013;Song 2016;Schmedes et al. 2013;Graves and Pitarka 2015;Crempien and Archuleta 2017;Withers et al. 2018).
We propose a kinematic earthquake rupture generator that combines the approach of Graves and Pitarka (GP2016) (Graves and Pitarka 2016) with the multiple asperity rupture characterization approach of Irikura and Miyake (IM2011) (Irikura and Miyake 2011). The resulting multi-scale hybrid rupture model incorporates distinct features of both original rupture models such as small-scale stochastic rupture variability and depth dependent rupture velocity and slip rate, inherited from GP2016, and distinct large slip To be submitted to special issue of Pure and Applied Geophysics, entitled ''Best Practices in Physics-based Fault Rupture Models for Seismic Hazard Assessment of Nuclear Installations''. patches required to produce observed near-fault local rupture directivity effects, coming from IM2011. The motivation of the work described here came from results of recent studies of near-fault ground motion from the M7 2016 Kumamoto, Japan earthquake Kawase et al. 2017;Kobayashi et al. 2017). The earthquake, which occurred in a densely populated area, caused significant damage in a relatively small area near the Futagawa fault (e.g. Kawase et al. 2017), one of the four fault segments activated during the earthquake. Several features of this earthquake were somewhat surprising and warrant further analysis. First, the extension of the earthquake damage zone does not correlate with the local site conditions (e.g. Kawase et al. 2017), suggesting that enhanced strong ground motion originated at the source. Second, in at least three near-fault sites the ground motion on the fault parallel component was as high as that on the fault normal component. This is inconsistent with the more typical case involving forward rupture directivity along a predominantly strike-slip fault rupture. Third, while successful at reproducing the time histories of acceleration and velocity, broadband simulations using rupture models generated with the IM2011 recipe underpredicted the permanent displacement at near-fault sites. The deficiency in long-period energy observed in simulated ground motion was caused by the absence of shallow slip (depth shallower than 3 km) for the entire fault, restricted by the recipe for generating the rupture model. Investigations of strong ground motion and rupture kinematics of the Kumamoto earthquake suggest that the rupture was dominated by a shallow area of large slip with relatively long rise time and a deeper area of relatively high slip rate . The longperiod energy radiated from the shallow rupture largely affected the near-fault permanent slip, and the high-slip rate area might have generated high frequency ground motion that was observed in the heavy damage zone (Kawase et al. 2017). These features are similar to rupture characteristics observed during recent large subduction zone earthquakes (e.g. Frankel 2017). Because large crustal earthquakes break the entire seismogenic zone, the rupture and, consequently, the radiated seismic energy is affected by depth-dependent variations of stress and frictional properties in the upper crust. This can explain the observed spatial separation between low-frequency and high-frequency generation areas, and variations of rupture kinematics with depth for large earthquakes. Such features affect the near-fault ground motion in a broad frequency range, and they need to be considered in refinements of kinematic rupture models for large crustal earthquakes.
In this study we investigated a modification of the original GP2016 rupture generator in order to produce slip models with stronger deterministic features, such as discrete areas with elevated slip. The modification allows for better control of the location and spatial extent of strong motion generation areas on the fault, while preserving original rupture features, such as depth dependency and correlation between rupture parameters, which are designed to reproduce expected effects of shallow slip on long-period ground motion. We analyzed the performance of the proposed method in broadband simulations (0-10 Hz) of near-fault ground motion from the Kumamoto earthquake, using the Graves and Pitarka (Graves and Pitarka 2010) simulation method. This was done by comparing simulated and recorded ground motion computed for three different rupture models, each representing distinct multi-scale features of the slip distribution for the Kumamoto earthquake. Comparisons with predicted ground motion using GMPEs were also used to investigate the effects of kinematic rupture model parameters on near-fault ground motion variability.
In the sections that follow, we first describe the approach to generate the hybrid model. We then describe the ground motion calculations for the Kumamoto earthquake that are used to examine the proposed methodology. These simulation results are also compared with estimates obtained from four NGA-West2 ground motions prediction equations (GMPEs). Finally, we show results from sensitivity tests of simulated near-fault ground motion to variations in the prescribed kinematic rupture parameters, such as rupture speed, slip duration rate, and slip contrast between the large slip areas and background slip.

Kinematic Rupture Models
Here we briefly describe the characteristics of the IM2011 and GP2016 rupture models and the process used to develop the proposed hybrid rupture model. IM2011 is based on the multiple-asperity concept of fault rupture and uses three sets of parameters, named outer, inner and extra fault parameters, to characterize the fault rupture kinematics. The outer parameters characterize the rupture area and magnitude, and the inner parameters define the spatial and temporal characteristics of slip distribution determined from estimated stress drop in the asperities and background areas of the fault. The extra fault parameters are the rupture nucleation location (hypocenter), rupture initiation point in each asperity, and rupture velocity. The outer and inner fault parameters are linked to the total seismic moment following empirical scaling laws. The number of asperities, total asperity area, and asperity slip contrast follows Somerville et al. (1999). In IM2011 the asperities are defined as rupture areas with both high slip (higher static stress drop) and high slip rate (faster rise time). This means that most of the strong shaking energy is generated in the asperity areas, which cover only a small part of the total fault area. In the case of the Kumamoto earthquake the asperity areas are confined to the deeper portions of the fault. They are separated from areas of largest slip located in the shallow part of the fault (e.g. Irikura and Kurahashi, 2018;Yoshida et al. 2017). This separation could be a typical feature of rupture kinematics for large crustal earthquakes (magnitude larger than 7). A similar trend in which high slip-rate areas are concentrated in the deeper part of the fault has also been observed during large subduction zone earthquakes (Kurahashi and Irikura 2013;Iwaki et al. 2016a, b). Kurahashi and Irikura (2013) refer to these high slip-rate patches as strong motion generation areas (SMGAs). Despite some similarities, the asperity areas adopted in the IM2011 should not be confused with deep subevents of subduction zone earthquakes that radiate most of the high frequency energy ([ 0.1 Hz) (e.g. Pulido et al. 2015;Frankel et al. 2018). The subevents do not contain the peak slip. This is the reason why they cannot be detected by long-period slip inversions. They are only detectable by broad-band modeling of ground motion waveforms (e.g. Frankel 2017;Wirth et al. 2017). More details on the IM2011 approach can be found in Irikura and Miyake (2011), Morikawa et al. (2011) and Pitarka et al. (2017).
The GP2016 rupture generator uses variable spatial and temporal kinematic rupture parameters that are calibrated using recorded ground motion and observed rupture kinematics. The rupture process, which is randomly heterogeneous at different scale lengths, controls coherent and incoherent interferences of waves generated at the source. The rupture generation process begins with the specification of a random slip field that is filtered to have a roughly wavenumber-squared falloff (e.g., Mai and Beroza 2002). The slip values are scaled to have a coefficient of variation of 0.85 and to also match the desired seismic moment. Given a prescribed hypocenter, the rupture propagation times across the fault are determined such that the average rupture speed scales at about 80% of the local shear wave velocity, and supershear rupture can happen in small areas of the fault. The rupture speed is further reduced by a factor of 0.6 for depths of 5 km and less, which is designed to represent the shallow, weak zone in surface-rupturing events (e.g., Marone and Scholz 1988;Dalguer et al. 2008;Pitarka et al. 2009). A perturbation is then applied to the rupture time at each subfault that is partially correlated with local slip such that the rupture tends to propagate faster in regions of large slip and slows down in regions of low slip. The slip-rate function is a Kostrov-like pulse (Liu et al. 2006) with a total duration (rise time) that is partially correlated with the square root of the local slip. Additionally, the rise time is scaled up by a factor of two within the 0-5 km depth range (Kagawa et al. 2004). The average rise time across the fault is constrained to scale in a self-similar manner with the seismic moment (Somerville et al. 1999). A complete description of the GP2016 approach can be found in Graves and Pitarka (2016).
The hybrid rupture generation process proposed here combines the discrete multiple asperity approach of IM2011 with the random spatial and temporal characterization of GP2016. The hybrid process begins with the generation of a slip distribution using GP2016. Then, one or more asperity patches are added to the GP2016 slip distribution. These asperities are rectangular in shape and the total number of asperities and their relative strength are guided by the rules given in the IM2011 approach. In specifying these asperities, the strength of each patch is given in terms of the average slip of the underlying GP2016 rupture. For example, an asperity with strength of one would have a slip equal to the average slip of the GP2016 rupture. Once the asperity patches are added to the GP2016 rupture, the slip values are rescaled to match the desired seismic moment. The remaining kinematic rupture parameters of rupture time and rise time are then determined using the original GP2016 process. The resulting hybrid rupture model contains both discrete high-slip rupture patches following IM2011 along with shorter lengthscale spatial and temporal heterogeneities following GP2016.
An example of this hybrid rupture generation process is shown in Fig. 1 for a hypothetical M 6.7 earthquake on a surface-rupturing vertical strike-slip fault. The initial GP2016 rupture is shown in the upper left set of panels of the figure. The subsequent panels of Fig. 1 show realizations of the hybrid approach where a single 10 km by 7 km asperity is added to the GP2016 rupture using scale factors of 1.0, 2.0 and 3.0, respectively. For each rupture, the figure displays both the slip and the peak slip-rate filtered at f \ 1 Hz. As the relative strength of the asperity is increased for the hybrid ruptures, more of the slip is concentrated in the asperity, and the level of background slip is decreased. Nonetheless, the effects of the randomized spatial and temporal heterogeneities from the GP2016 rupture are still apparent in the resulting hybrid rupture model. Additionally, not only does the addition of the asperity affect the slip distribution it also has an impact on the peak slip-rate of the rupture, with the level of the peak slip-rate increasing as the strength of the asperity is increased.

Simulation of M7 2016, Kumamoto, Japan Earthquake
The Kumamoto earthquake is a strike slip event that occurred on four fault segments (e.g. Yoshida et al. 2017). Figure 2 shows the location of the fault segments and strong motion stations used in our simulations. Three fault segments ruptured the surface, while the fourth, at the eastern end of the rupture, did not. The fault geometry and location are based on the distribution of aftershocks that occurred immediately after the main event, and InSAR (Geospatial Information Authority of Japan 2016) data. The orientation of the fourth, and most eastern fault segment, is not well resolved. Additionally, as suggested by fault slip inversions of long period waves, the amount of slip in this segment and its contribution to overall ground motion is relatively low ).

Rupture Models
The rupture models of the Kumamoto earthquake, tested here, are shown in Fig. 3a. The models share the same fault geometry, mechanism, and hypocenter location obtained from kinematic rupture inversions (e.g. Yoshida et al. 2017). The fault parameters are listed in Table 3. The first model, named IMh1 contains a fully deterministic slip, dominated by two distinct asperity areas with high slip. The size and depth of the asperity areas, and the contrast between the asperity slip and background slip are determined following the IM2011 recipe for a M7 strike slip fault. The location of the asperities was based on kinematic rupture models proposed for the earthquake Yoshida et al. 2017). In order to be consistent with the GP2016 model, the rupture speed for IMh1 is set at 80% of the local shear wave velocity, which is slightly higher than the value of 72% prescribed by the recipe. Additionally, there are no stochastic variations in slip, rupture time or rise time, consistent with the IM2011 approach. In the second model, named GP, the slip is fully stochastic. GP was selected among several rupture scenarios generated using the proposed method. The main requirement in the selection criteria was that, similar to the observed slip ), each fault segment should contain an area of elevated slip. The third rupture model, named HB2, is a combination of IMh1 and GP models, generated with the proposed hybrid procedure. Note that, as in the GP2016 method, our rupture generator provides partial correlation between rise time and the square root of local slip. This results in a tendency for the rise time to lengthen and the rupture speed to increase as the slip increases. Also apparent is the systematic reduction of rupture speed and lengthening of rise time along both the top (upper 5 km) and bottom portions of the rupture, as dictated by the GP2016 Figure 1 Example rupture models generated using the GP2016 approach (upper left), and three realizations using the hybrid approach: HB1.0 (upper right) adds a slip patch scaled at 1.0 times the average slip of GP2016, HB2.0 (bottom left) adds a slip path scaled at 2.0 times the average GP2016 slip, and HB3.0 (bottom right) adds a slip patch scaled at 3.0 times the average GP2016 slip. For each rupture, the top panel displays the slip distribution along with rupture time contours at 2 s intervals, and the bottom panel displays the distribution of peak slip-rate filtered at f \ 1 Hz. The triplet of numbers above each panel gives the minimum, average and maximum value, respectively, of the parameter shown in the panel Vol. 177, (2020) Kinematic Rupture Modeling of Ground Motion from the M7 Kumamoto, Japan Earthquake 2203 approach. The GP and HB2 models shown in Fig. 3a were generated using default values for the kinematic rupture parameters, as described in Graves and Pitarka (2016). These distinct features are consistent with observed rupture characteristics of the Kumamoto earthquake. For example, several studies have shown that the Kumamoto earthquake rupture had at least one shallow area of large slip, and two deeper areas of large slip rate (e.g. Yoshida et al. 2017;Asano and Iwata 2016;Irikura and Kurahashi 2018). A similar trend in which high slip rate areas are concentrated in the deeper part of the fault has also been observed during large subduction zone earthquakes (Frankel 2017;Kurahashi and Irikura 2013). These features for ruptures that break the entire seismogenic zone, can be explained by differences in fault rupture dynamics caused by difference in frictional properties between materials in the shallow part (weak zone) and deeper part of the crust (e.g. Pitarka et al. 2009;Dalguer et al. 2008;Song 2016;Dan et al. 2018). The depth scaling of rise time and rupture speed prescribed in the GP2016 approach also reflects key observations where large shallow fault slip does not necessarily translate into large high-frequency motion radiation. This means that strong radiation of shorter period motion does not necessarily coincide with regions of large slip. These two important features are embedded in the HB2 rupture model.

Ground Motion Modelling
The ground motion simulations were performed using the Graves and Pitarka hybrid simulation procedure (Graves and Pitarka 2010). We compute free-surface broadband (0-10 Hz) ground motion at 19 stations surrounding the fault, extending to a fault distance of about 32 km. Table 1 lists the station locations and local V s30 . The V s30 at each station was based on personal communications with Somei (2018). We used the three rupture models discussed in the previous section to run the simulations. The sub-fault dimensions used in computing the low frequency part of ground motion were 0.2 km by 0.2 km, and 2 km by 2 km in computing the high frequency part of ground motion. The fault depth was set at 0 km. In the simulations of the high frequency part of ground motion we used the default stress parameter of 50 bars, adopted in simulations of earthquakes in California. The low-frequency part of ground motion (0-1 Hz) was calculated using synthetic Green's functions computed with the FK method of Zhu and Rivera (2002) for a prescribed 1D velocity model, listed in Table 2. The velocity model was obtained by averaging the 1D velocity models extracted from the NIED 3D velocity model (Fujiwara et al. 2012) at each of the 19 stations considered in this study. The deeper structure used beneath seismic bedrock (V s = 3.4 km/s) in the NIED model was modified from the 1D model of Wu et al. (2008). The V s 30 of the prescribed 1D structure is set at 500 m/s. The broadband synthetic seismograms were corrected for site effects using local V s 30 and site response adjustment factors taken from the GMPE of Campbell and Bozorgnia (2014) and implemented using the approach of Graves and Pitarka (2010).
A synthesis of the simulations' performance is shown in Fig. 3b. We show the bias of RotD50 The residuals used to determine the goodness of fit are computed between the simulations and the recorded data Vol. 177, (2020) Kinematic Rupture Modeling of Ground Motion from the M7 Kumamoto, Japan Earthquake (Boore 2010) pseudo spectral acceleration response with 5% damping (ln(Rec/Syn)) between recorded (Rec) and synthetic (Syn) ground motion, averaged over 19 stations, also known as goodness-of-fit plot (Goulet et al. 2015). Overall the bias obtained with the three models is small. This result demonstrates that the rupture models match the general characteristics of the recorded near-fault motion equally well, in a broad frequency range (0-10 Hz).
IMh1 performs well in reproducing the observed spectral response. For this model we used the GP rupture model default values for both rupture speed V r and stress parameter. V r is set at 80% of local shear-wave speed V s , and the stress parameter is set at 50 bars. In ground motion simulations these two parameters trade off with each other at high frequencies (f [ 1 Hz). For example, if in generating the IMh1 model, we use V r = 0.72*V s , as recommended Table 1 List and location of stations  Comparison of free-surface recorded (black traces) and simulated (red traces) broadband (0-10 Hz) ground motion time histories of acceleration, velocity and displacement computed with the HB2 rupture model shown in Fig. 3 Vol. 177, (2020) Kinematic Rupture Modeling of Ground Motion from the M7 Kumamoto, Japan Earthquake by the IM2011 recipe for this type of rupture and magnitude, a higher stress parameter of 60 bars is required to produce an equivalent fit to the spectral response. The difference in stress parameter values between the GP and IM2011 models is due to their different technical approaches followed in computing the high frequency part of ground motion. For large earthquakes the stress parameter recommended by IM2011 is as low as 31 bars (Fujii and Matsu'ura 2000).
Although not analyzed in detail here, another factor that trades off with rupture velocity in affecting the computed ground motion is the Vs30-based site response correction. It is possible that the empirical site correction model of CB2014 that is used here may not work well for site corrections in Japan. Consequently, at this point we hesitate to draw a firm conclusion about the increased stress parameter required in simulations of the Kumamoto earthquake when using the IM2011 recipe based rupture velocity Vr = 0.72*V s .
The GP model, which is mainly calibrated for California earthquakes, slightly under predicts the ground motion in the short-period range of 0.1-1 s. This small discrepancy is eliminated in the simulation with the HB2 model. By design, the HB2 model better characterizes deep large slip rate areas, affecting short period seismic energy (\ 1 s), and shallow slip areas, affecting periods longer than 1 s. Note that, although the bias for both HB2 and IMh1 is very low for all response periods, the standard deviation obtained with HB2 is higher for periods longer than 1 s. This indicates that compared to GP the proposed rupture model may produce relatively higher intraevent ground motion variability. This is a desired feature since in general broadband simulation methods underestimate the intra-event ground motion variability (e.g. Crempien and Archuleta 2017).
Figures 4, 5 and 6 illustrate the comparison of recorded and simulated free-surface acceleration, velocity and displacement time histories, computed at 19 stations using the IMh1, GP, and HB2 rupture models, respectively. As already indicated by the goodness-of-fit plots, with some exceptions, the waveforms fit obtained with the three models is acceptable even at relatively long distances. Note that at some stations (e.g. KMMH09, KMM006) the recorded displacement has a late-time drift caused by the instrument response. We focused our attention on the simulated waveforms of velocity and displacement at two near-fault stations, KMMH16 and KMM005, located 2 km and 5 km off the fault trace, respectively. KMMH16 located within the heavy damage zone, and recorded the highest ground motion acceleration and velocity (Kawase et al. 2017). The waveform comparison between the recorded and synthetic ground motion is shown in Fig. 7. The effects of the IMh1 model on near-fault ground motion are mainly controlled by rupture directivity caused by the uniform rupture of the two large asperity areas with constant high slip. These effects are manifest by increased ground velocity in the propagation direction of the rupture in a narrow zone along the fault. As it will be shown in a subsequent section, the planar fault assumption and fault dipping toward north enhance rupture directivity effects in an area north of the fault trace where both KMM005 and KMMH16 are located. KMMH16 appears to be located near the edge of the amplification zone. The local rupture directivity effect is sensitive to the relative location of the station with respect to asperity and hypocenter. For these reasons IMh1 reproduces the observed peak ground motion velocity (PGV) at KMM005, but underestimates the PGV and peak displacement at KMMH16.
We concluded that IMh1reproduces the general characteristics of the recorded ground motion reasonably well. IMh1 has the potential to generate strong ground motion with directivity effects at areas near the asperities. The spatial extent and location of the amplification zones caused by rupture directivity effects depends on the asperity and rupture starting point relative locations (e.g., Somei 2018). However, due to reduced amount of near-surface slip, enforced by the IM2011 recipe, IMh1 underestimates long period ground motion and permanent displacement observed at a near-fault sites. The GP rupture model does a good job at reproducing the ground motion amplitude at KMM005 and KMMH16, including permanent displacement on the N-S component. At these two nearfault sites a similar performance was obtained with the hybrid model HB2. Both models reproduce the recorded velocity reasonably well. Their performance is slightly better than that of IMh1 in the sense that at station KMMH16 they generate larger ground motion velocity, and much larger permanent displacement than IMh1. We concluded that the lack of an extended area with near-surface slip, not well represented by our slip models, and potential local changes of the focal mechanisms, not included in our models, caused the mismatch of the large permanent displacement on the E-W and vertical components at KMM005 and KMMH16 stations (e.g., Somei 2018) ( Table 3).

Sensitivity of Ground Motions to Rupture Parameters in HB2 Model
We investigated the effects of HB2 rupture parameters on simulated near-fault ground motion   modified values of rupture parameters are shown in Table 4. Figure 8 shows results of these simulations displayed in two ways. The first uses the goodness-of-fit plots, averaged over 19 recording sites, that we have already discussed. The second uses maps of faultnormal (FN) and fault-parallel (FP) of peak ground motion velocity, and their ratio, computed for each rupture scenario on a dense grid of 528 stations, covering the entire area shown in Fig. 2. The stations grid spacing is 4 km.
Increasing V r /V s from 0.7 (recommended in IM2011) to 0.85, which is slightly above the value of 0.82, found for the Kumamoto earthquake , results in increased ground motion amplitude on a broad period range. The corresponding ground motion maps indicate that a faster rupture can double the area of strong ground motion. The location of the area of largest ground motion velocity, located between stations KMM005 and KMMH16, remains unchanged. Based on the goodness-of fit plots we conclude that an average V r /V s = 0.7 systematically under predicts the level of observed ground motions and is, therefore, too low for this earthquake. The sensitivity of ground motion to rupture velocity and its spatial variation has been the subject of several studies (e.g., Iwaki et al. 2016a, b).
The effect of increasing the average rise time duration, which in our rupture model is controlled by the risetime coefficient (Graves and Pitarka, 2015), is significant only in a narrow period range, between 1 and 2 s. For those periods the increase in rise time coefficient from 0.8 to 3.2 reduces the ground motion amplitude by about 50%. The peak amplitude reduction can be explained by the fact that for the same seismic moment, the increase of rise time decreases the peak slip-rate, which in turns reduces the amplitude of ground motion velocity and acceleration. In our analysis the rupture scenario with longer rise time produces a smaller area of strong ground motion. The goodness-of-fit plots obtained for risetime coefficients of 0.8 and 3.2 suggest that these values, which are below and above the default value of 1.6, respectively, are acceptable, and that variations of risetime coefficient bounded by these values only affect ground motion at periods between 1 and 2 s.
The third rupture parameter considered in these analysis is the asperity slip ratio, measured as the ratio between the average slip in the asperities and background areas. The effect of increasing the asperity slip ratio from 2, recommended by the IM2011 recipe, to 4, on the ground motion goodnessof-fit and spatial distribution of PGV are shown in Sensitivity of spectral acceleration goodness-of-fit plot and PGV maps to asperity slip ratio in simulations of the 2016 Kumamoto earthquake using a realization of the HB2 rupture model with an asperity slip ratio of 4 (top panel) and default slip asperity ratio of 2 (bottom panel). Note that in the simulations shown in Fig. 8 we used the default value of 2 for the asperity slip ratio Vol. 177, (2020) Kinematic Rupture Modeling of Ground Motion from the M7 Kumamoto, Japan Earthquake Fig. 9. The ground motion amplitude is reduced by at least 30% in the long period range of 3-7 s. This can be explained by the reduction of background slip, and in particular, the shallow slip, required to compensate for the slip increase in the asperity areas, while preserving the target seismic moment. Because this is a long-period effect the ground motion amplitude reduction is not observed on the PGV maps.
In the final stage of our study we used the ground motion predicted by the GMPEs as a reference in comparing effects of rupture parameters, described above, on simulated intra-event ground motion variability. First, we compare the GMPEs with recorded and simulated ground motion for the HB2 rupture model, as shown in Fig. 10a, b. The GMPEs were computed for a Vs30 = 300 m/s and the synthetic ground motions are all corrected for site effects using local V s30 and site response adjustment factors, taken from the GMPE of Campbell and Bozorgnia (2014). Except for 5 s period for which the recorded ground motion at distant stations is smaller than the one predicted by the GMPEs, the recorded RotD50 PSA is similar to the average value of the GMPEs. In addition, the simulated ground motion fits the GMPEs equally well at all periods. These favorable comparisons suggest that the near-fault ground motion generated during the Kumamoto earthquake is similar to the average ground motion predicted by the GMPEs for a strike slip M7 earthquake. Therefore the GMPE predictions can be used as a reference in analysis of effects of HB2 rupture parameterization on computed ground motion. Figure 11 shows the comparison of ground motion predicted by the GMPEs and the simulated ones using V r /V s ratios of 0.7, 0.8, and 0.85, respectively, where 0.8 is the default value. The effects of rupture speed variation are significant only for spectral periods of 0.2 s and 0.3 s, which correspond to frequencies of 5 Hz and 3 Hz. At these frequencies increasing the rupture speed from Vr = 0.7*Vs to Vr = 0.85*Vs amplifies the ground motion by as much as a factor of 2. The effects are present, but less pronounced at longer periods.  Figure 12 shows the comparison of ground motion predicted by the GMPEs and the simulated ones using risetime coefficients of 0.8, 1.6, and 3.2, respectively. As mentioned above, the rise time is a parameter of the slip rate function, which in our broadband simulation method is only used in computing the long-period part ([ 1 s) of ground motion. Therefore, effects of changing the rise time are only Figure 11 Comparison of ground motion predicted by the GMPEs and the simulated ones (gold circles) using HB2 rupture model for V r /V s ratios of 0.7, 0.8, and 0.85, respectively. V r /V s ratio is indicated on top of each panel Vol. 177, (2020) Kinematic Rupture Modeling of Ground Motion from the M7 Kumamoto, Japan Earthquake 2215 expected for periods longer than 1 s. We observe that at periods of 2 s and 5 s the simulated spectral response is systematically reduced by 20% and 100%, as the risetime coefficient increases from 0.8 to 1.6, and 3.2, respectively. This result suggests that a risetime coefficient of 3.2 may generate ground motions that are too low for periods longer than 2 s. We also looked at the impacts of rupture speed and rise time on intra-event ground motion variability. This is done by computing the standard deviation of simulated ground motion (sigma) with respect to the average GMPE, and analyzing its variability as a function of period and distance. We grouped the stations into four distance bins so that we could perform the statistical analysis of sigma. The distance bins are 2 km, 5 km, 10 km, and 20 km, and number of stations in each bin are 8, 42, 50, and 67, respectively. For each rupture realization, we computed the standard deviation (sigma) of simulated spectral acceleration across all stations within each distance bin at each oscillator period. This is a measure of intra-event variability since it considers the motions for only a single rupture realization. Figure 13 shows the simulated intra-event sigma computed for three different rupture speeds, Vr = 0.7*Vs, Vr = 0.8*Vs, and Vr = 0.85*Vs. A zero bias means that the ground motion averaged over the stations is the same as the average ground motion predicted by the GMPEs.
In terms of ground motion amplitudes, the main conclusion derived from the sigma plots is that the ground motion increases with increasing rupture speed at all spectral response periods. In terms of ground motion variability, we found that not only does it depend on the spectral period, but it also depends on distance from the fault. This indicates that wave propagation scattering impacts the frequency content of the signals as the waves propagate away from the fault. The general trend seen in the simulated ground motion variability is that at near-fault distances (distance = 2 km) the scenario with the lowest rupture speed produces strong period dependent variability. At distances longer than 2 km the variability becomes more uniform, and is similar among the three rupture speeds. It is important to note that, as indicated by the shadowed yellow zone, the standard deviation among the four GMPEs greatly increases with period at all fault distances. A similar trend is observed in the synthetic data for the three rupture speeds. The impact of rise time on intra-event ground motion variability is shown in Fig. 14. As with rupture speed, we computed the standard deviation of simulated ground motion (sigma) with respect to the average GMPE, and analyzed its variability for three rupture scenarios in which the risetime coefficient is 0.8, 1.6, and 3.2, respectively. The impact of slip-rate function duration on the ground motion variability is smaller than that observed for the rupture velocity.
Our most significant observation is that the variability is largest at the fault distance of 2 km. At this distance the ground motion variability has the tendency to increase with period. This is in agreement with the trend seen in the standard deviation among the four GMPEs (shaded area in the sigma plots). The variability becomes less period dependent at longer distances for the three rise times considered here.

Conclusions
In this study we examined the performance of a kinematic rupture generator that combines the randomized spatial field approach of Graves and Pitarka (2016) with the multiple asperity characterization approach of Irikura and Miyake (2011)  simulations of the M7, 2016 Kumamoto, Japan earthquake. The resulting hybrid rupture model contains both discrete high-slip patches following IM2011, along with shorter length-scale spatial and temporal heterogeneities coming from GP2016. Except for allowing slip concentration in rectangular asperities, the proposed model has the same features as the GP2016 model from which it was derived. Those features are the depth dependency of rupture speed and rise time, and correlated small scale stochastic variability of all rupture parameters. We used three rupture models, named IMh1, GP, and HB2, each representing distinct rupture features produced by the rupture generator, to simulate ground motion observed at 19 stations during the M7, 2016 Kumamoto Japan earthquake. The simulations were performed using the broadband simulation technique of Graves and Pitarka (2015). The overall performance of the three rupture models at reproducing characteristics of the recorded motion across the considered frequency range of 0-10 Hz is generally acceptable. The IMh1 model, in which the fully deterministic slip distribution follows the IM2011 recipe, tends to underestimate the long period ground motion and permanent displacement at a near-fault sites. This discrepancy is caused by the absence of the near-surface slip enforced by the recipe. Recent proposed improvements of the recipe ) allow for inclusion of areas of shallow slip in addition to strong motion generation areas.
GP, a fully stochastic model, slightly under predicts the ground motion in the short-period range of 0.1-1 s. This small discrepancy is eliminated in the simulation with the hybrid model, referred as HB2 model. Designed to better capture the deterministic features of the slip distribution, such as discrete areas with large slip rate, and shallow slip areas with longer rise time, HB2 generates larger long period ground motion than IMh1. We recommend that in order to correctly produce the large near-fault permanent displacement generated by large crustal earthquakes, similar to the Kumamoto earthquake, rupture models should contain a large slip zone that reaches the free surface. This feature can be accommodated by the proposed hybrid generator which was designed to better control the location and size of areas with large slip.
We also investigated the sensitivity of simulated ground motion to parameters of the HB2-type rupture scenarios of the Kumamoto earthquake. The impact of increasing the average rupture speed ratio (V r /V s ) from 0.7 (recommended in IM2011) to 0.85, which is slightly above the value of 0.82 found for the Kumamoto earthquake, results in increased and ground motion amplitude on a broad period range. Maps of simulated PGV indicate that a faster rupture can double the area of larger ground motion. We concluded that an average V r /V s = 0.7 is too low for the Kumamoto earthquake. The effect of increasing the average rise time duration is significant only in a narrow period range, between 1-2 s. For those periods doubling of rise time duration reduces the ground motion amplitude by 50%. We also found that increasing the slip amount in asperity areas results in reduction of long period ([ 5 s) ground motion amplitude, primarily due to the reduction in average slip in the background areas.
In our sensitivity analysis of intra-event ground motion variability to rupture speed and rise time, we found that near-fault ground motion variability is period dependent. In general, at very near-fault distances (distance = 2 km) the scenario with the lowest rupture speed produces the strongest period dependent variability. The impact of rise time on ground motion variability is smaller than that observed for the rupture velocity. Our simulations indicate that at 2 km fault distance the ground motion variability has a slight tendency to increase with period. At distances 5 km and longer, the variability becomes less period dependent for any of the rupture velocities and rise times considered here.
the Nuclear Regulation Authority (NRA), Japan.

LLNL-JRNL-764062
Open Access This 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.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.