Seismological asperities from the point of view of dynamic rupture modeling: the 2007 Mw6.6 Chuetsu-Oki, Japan, earthquake

We study the ground motion simulations based on three finite-source models for the 2007 Mw6.6 Niigata Chuetsu-oki, Japan, earthquake in order to discuss the performance of the input ground motion estimations for the near-field seismic hazard analysis. The three models include a kinematic source inverted from the regional accelerations, a dynamic source on a planar fault with three asperities inferred from the very-near-field ground motion particle motions, and another dynamic source model with conjugate fault segments. The ground motions are calculated for an available 3D geological model using a finite-difference method. For the comparison, we apply a goodness-of-fit score to the ground motion parameters at different stations, including the nearest one that is almost directly above the ruptured fault segments. The dynamic rupture models show good performance. We find that seismologically inferred earthquake asperities on a single fault plane can be expressed with two conjugate segments. The rupture transfer from one segment to another can generate a significant radiation; this could be interpreted as an asperity projected onto a single fault plane. This example illustrates the importance of the fault geometry that has to be taken into account when estimating the very-near-field ground motion.


Introduction
Predicting the ground motions for any given scenario of an earthquake is an important seismological task for seismic hazard evaluation (Douglas and Aochi 2008). Nowadays, the earthquake model is kinematically constructed based on the statistical analyses of past earthquakes obtained from various inversions or synthetic earthquake scenarios dynamically simulated (e.g., Mai and Beroza 2002;Irikura and Miyake 2011;Song et al. 2013). In some research projects, there have been attempts to carry out the ground motion simulations using the dynamically simulated earthquake models directly (e.g. Olsen et al. 2009). Indeed, for two decades, dynamic rupture models have been more commonly applied to reproduce the ground motions of recent earthquakes, such as the 1992 Landers earthquake (Olsen et al. 1997;Peyrat et al. 2001;Aochi et al. 2003). Furthermore, the characteristics of the ground motion based on the dynamic rupture models are synthetically studied in terms of rupture velocity, fault geometry, heterogeneity, and so on (e.g., Oglesby and Day 2002;Aochi and Olsen 2004;Aochi and Douglas 2006;Schmedes and Archuleta 2008;Dunham and Bhat 2008). Regardless of the progress in dynamic rupture modeling, the models are difficult to produce and calibrate. In this paper, we aim to show the applicability and the utility of the dynamic rupture models for the nearfield ground motion simulations applied to the 2007 Mw6.6 Chuetsu-oki, Japan, earthquake.
The 16th July 2007 Chuetsu-oki earthquake (Niigata prefecture, Japan; Fig. 1) led to some damage in Kashiwazaki and the shutdown of the Kashiwazaki-Kariwa Nuclear Power Plant (NPP), located near station KSH. The dense seismological observational networks, the geodetic observations, and the aftershock analyses in this region reveal the complexity of this earthquake in every aspect: source, propagation, and site (e.g., International Atomic Energy Agency 2007). The capability of reproducing the observed records through numerical modeling becomes a key issue for testing our knowledge and understanding of the quantitative seismic hazard assessment at any location, in particular in the very near field of the causal earthquake source area (e.g., Irikura 2010). This paper focuses on the input ground motion at the nearest seismological observations-the borehole data within the site of Kashiwazaki-Kariwa NPP (KSH-SG4). Utilizing the best available 3D geological model, we carry out the ground motion simulations based on different earthquake source models including both kinematic and dynamic descriptions. The purpose was to emphasize how kinematically inferred earthquake asperities look from the dynamic rupture models and to discuss the applicability of the dynamically simulated source models for estimating input ground motion.  Sekiguchi et al. (2009). Triangles represent the permanent seismological observation networks, and the names are given for those which are going to be used in this study. The acceleration records are at the two nearest stations, KSH (SG4 at depth of 250 m) and NIG018. No frequency filter is applied 2 Source models of the 2007 Chuetsu-Oki earthquake Figure 2 summarizes three source models for the 2007 Chuetsu-oki earthquake used in this study. Various researchers used inversion to obtain kinematic descriptions of the source process. The fault orientation (SEor NW-dipping thrust faulting) was not evident just after the earthquake mainly because the estimated rupture zone is off the coast and the routinely determined location of the aftershocks were inaccurate (Aoi et al. 2008). The observed strong ground motion along the coastline suggested a NW-dipping fault such that the ground motion could be influenced by rupture directivity. However, the interpretation for faulting based on a SEdipping faulting became widely accepted (e.g., Miyake et al. 2010) according to the relocation of the aftershocks from temporary observation networks (e.g., Shinohara et al. 2008;Kato et al. 2008). For this study, we adopt one of the SE-dipping fault models, model B of Aoi et al. (2008), as a reference, though we refer to it as model K (Fig. 2). This model was also used in previous ground motion simulations (Aochi et al. 2013a).
The strong ground motion observed in the near field (a distance from the ruptured fault plane closer than 10 km) was the subject of a major debate about the  Aoi et al. (2008). Model AD: a dynamic source model with three strong motion generation areas simulated in Aochi and Dupros (2011). Model AK: a dynamic source model with rupture on two cross-cutting faults (NW-and SE-dipping segments in the north and south, respectively) simulated in Aochi and Kato (2010). Stars represent the epicenter locations in the map. In the bottom, the rupture time (T RUP ), maximum slip rate, and final slip distribution are shown, respectively, with the same scale. For the points of the maximum slip in the three models (denoted by a cross), slip rate evolutions are illustrated at the bottom reliability of seismic hazard estimation (Strasser and Bommer 2009;Baumann and Dalguer 2014). The nearest seismic stations (KSH and NIG018;Figs. 1 and 2) are located above the inferred ruptured area along the fault strike. The ground motion at these stations is significantly affected by near-field effects of the earthquake source, such as asperity locations. They were not used in the inversion of Aoi et al. (2008). Irikura (2008) proposed a characteristic model consisting of three asperities, identified as strong motion generation areas (SMGAs). These localized asperities should produce recognizable phases in the seismograms. This characteristic earthquake model was further studied by dynamic modeling assuming a slip-weakening friction law (Aochi and Dupros 2011) and re-illustrated in Fig. 3. From the dynamic point of view, these SMGAs correspond to areas that have stress drops two times larger than the rest of the fault area. This is consistent with the generalized strong motion recipe of Irikura and Miyake (2011). Irikura (2008) proposed that the rupture of the third asperity located in the southernmost end was in the opposite direction of the general southward rupture propagation of this earthquake (multi-hypocenter model). In order to reproduce this effect dynamically, the third asperity needs to be partially surrounded by barriers so that the rupture first moves to the southern edge of the asperity and then ruptures to the northern edge (see the snapshot at 8.3 s in Fig. 3). This dynamic model is presented as model AD.
Most geodetic studies, on the other hand, propose two segmented, conjugate fault planes (Aoki et al. 2008;Nishimura et al. 2008), which were also inferred from the most accurate aftershock distributions (Shinohara et al. 2008;Kato et al. 2008). Aochi and Kato (2010) previously carried out dynamic rupture simulations for the conjugate fault geometry consisting of a northern NW-dipping segment (segment N) and a southern SEdipping one (segment S). They varied the intersection angle and the overlapping distance, as illustrated in Fig. 4. Among their 45 simulations, only 14 show a sequential rupture transfer from one to another segment. We adopt one scenario in which the model parameters are moderate: a frictional coefficient of 0.3, an overlapping distance of two segments of 4 km, and a dipping angle of SE-dipping fault of 45°. This is designated as model AK.
3 Numerical simulations of ground motion

Simulation methods
The complexity of the 3D geological structure is well known in this region due to its complex tectonics. Aochi et al. (2013a) used the three available structure models of this region Fujiwara et al. 2009;Sekiguchi et al. 2009) for ground motion simulations and discussed the validity of each model. Ground motion using the 3D structure obtained by tomography ) shows a good coherence on rock sites that are situated close to each other. On the other hand, the models calibrated based on the geological map and the geophysical cross-sections (Fujiwara et al. 2009;Sekiguchi et al. 2009) are generally suitable for the soft sites. However, even when the models are good enough for the direct waves at certain stations, it is still difficult Fig. 3 Dynamic model parameters (initial shear stress, τ 0 ; peak strength, τ p ; dynamic stress, τ r ; and critical slip displacement, D c ) and snapshots of slip velocity and shear stress for source model AD-dynamically simulated as rupture on a single fault plane with three asperities (modified after Aochi and Dupros 2011). Note that the slip-weakening rate, (τ p − τ r )/D c , is the same everywhere to obtain coherent phases for the later seismic arrivals. In this study, we adopt the 3D structure model by Sekiguchi et al. (2009) which has integrated geological layers, as shown in Fig. 1, and has a minimum shear velocity of V min = 400 m/s. For ground motion simulations, we use a fourth-order finite-difference method (Aochi et al. 2013a, b) with a grid spacing of Δs = 80 m for a dimension of 110 km (EW) × 120 km (NS) × 30 km (UD). The ground surface is approximated as flat, and the Sea of Japan is not taken into account. The maximum reliable frequency is estimated as f max = V min /(5 ⋅ Δs) = 400/(5 ⋅ 80) = 1.0 Hz. The time step is 0.004 s; the calculation is run for a duration of 60 s. The simulation procedure for generating ground motion is basically the same for both the kinematic and dynamic source models, as in Aochi and Dupros (2011) and Aochi et al. (2013a), but using finer grids in this study. Any finite source can be introduced as a series of point sources with a predefined slip function of any arbitrary shape (Aochi et al. 2013b).
The included source model has been previously presented; its characteristics are summarized in Fig. 2. We propose the following models: the kinematic model K has a large rupture area, while the dynamic rupture models AD and AK have shorter rupture dimensions. The slip velocity function is also smooth (lower peak and longer duration) in model K and is shaper in AD and AK.

Simulation results
Figures 5, 6, and 7 show the simulated ground motions for the nearby stations. For model K, the synthetic ground motions have been presented by Aochi et al. (2013a). The dynamic rupture models (models AD and AK) were computed using the 3D geological model. The ground motions are aligned at the origin time of each simulation which coincides with the origin time of the hypocenter (10:13:22.16 local time). All ground motions are bandpass-filtered between 0.1 and 1.0 Hz. The kinematic model K uses a simpler 1D structure. Using only a portion (14 s around the direct S wave's arrival) of the seismograms, the reconstruction of the full waves in a 3D structure does not always ensure the coherent fitting of the waveforms even within this frequency range. At the nearest station, NIG018, it is known that the soft soil had been subjected to Fig. 4 Dynamic model parameters (initial shear stress, τ 0 ; peak strength, τ p ; dynamic stress, τ r ; and initial normal stress, σ 0 . Note that the slip-weakening rate is 30 MPa/m and is given everywhere) and snapshots of the slip rate and shear stress for the dynamic rupture process (model AK) on two conjugate fault segments (modified after Aochi and Kato 2010) liquefaction ( Fig. 1) so that the synthetic ground motions cannot be compared directly to the observations.
It is difficult to capture the waveform characteristics precisely in the backward direction of the rupture propagation, namely, in the northeast area (NIG10 and NIG11). However, in the forward direction in the southwest area (NIG024 and NIG025), we find that the characteristic waveforms are better simulated by the dynamic model (AK) than by the kinematic model (K). This is an important feature. The kinematic model is based on an inversion that might have missed or underestimated significant behavior of the rupture process due to a priori constraints or smoothing of the inversion. Therefore, the source model might underestimate the strong ground motions that are closely related to the source process. The dynamic models AD and AK produce similar waveforms, in particular at NIG024, NIG025, NIG004, and NIG019, and seem to be closer to the observations than model K. The similarity of the resultant ground motions for the two models, AD and AK, indicates that the Comparison of the synthetic NS ground motion from the three source models (K, AD, and AK) with the recorded NS velocity. The ground motions were filtered between 0.1 and 1.0 Hz. For reference, the source model K is illustrated on the map asperities inferred on a single fault plane can be represented by the geometrical irregularities of the fault system.
In Fig. 8, we compare two additional simulations in which we use either segment (N only or S only) of model AK to clarify the contribution of the non-planar fault structure. At the closest KSH station, the contribution is clear. The main features of the waveforms are reproduced by segment S, which is closer to KSH and releases 71 % of the total seismic energy. On the other Fig. 7 Comparison of the synthetic UD ground motion from the three source models (K, AD, and AK) with the recorded UD velocity. The ground motions were filtered between 0.1 and 1.0 Hz. For reference, the source model K is illustrated on the map Fig. 8 Comparison of the synthetic ground motions from either segment (segment N or S only) of model AK with the recorded velocity. The ground motions were filtered between 0.1 and 1.0 Hz. For reference, model AK (namely, the summation of segments N and S) is also shown hand, it is observed that the first pulse is represented by segment N. Such identification becomes difficult at NIG004 as the wave propagation is perturbed during a longer propagation distance. We also compare the final displacement field on the ground surface among models K, AD, and AK in Fig. 9 since two fault segments are proposed by geodetic observation (Aoki et al. 2008;Nishimura et al. 2008). A difference is hardly found in the vertical component (UD). In the horizontal components, model AK generates a clearer change in trend, although, again, the field is governed principally by segment S. The comparison between the planar models (K and AD) and the non-planar one (AK) has important implications for seismic hazard. Given that it is difficult to identify a priori any asperities on a fault plane before an earthquake, it may be more important to determine the fault structure first.

Criteria for ground motion estimation
This paper does not aim to improve the model parameters, but tries to show the performance of different earthquake models in generating the ground motion, in particular the performance of the two dynamic rupture models. Although inversion of dynamic rupture parameters has been done for a decade (Peyrat et al. 2004;Ruiz and Madariaga 2011;Douilly et al. 2015), the number of inverted model parameters is usually limited to about 10 due to the high nonlinearity of the system and computational costs. This means that one cannot expect the same degree of spatial or temporal resolution as found in kinematic inversions unless the model parameters of dynamic rupture are constrained by kinematic inversion results (e.g., Peyrat et al. 2004). For the purpose of the ground motion prediction for seismic hazard, the engineers are more interested in the ground motion parameters than in coherent waveforms. Aochi and Douglas (2006) proposed to compare statistically the ground motion parameters at numerous points from the simulations with the ground motion prediction equations in terms of peak ground acceleration (PGA), peak ground velocity (PGV), response spectral acceleration, Arias intensity, and relative significant duration. Olsen and Mayhew (2010) propose a goodness-of-fit (GOF) criterion for the validation of the broadband synthetics, which consists of different tests on ground motion metrics. The GOF score (normalized to 100) is defined as the complementary error function (erfc) of a normalized residual, NR (Olsen and Mayhew 2010) By associating a weight on the GOF score calculated for each of the metrics, the average GOF score is obtained. Olsen and Mayhew (2010) used the metrics consisting of PGA, PGV, peak ground displacement (PGD), averaged response spectral acceleration (RS), Fourier spectrum (FS), energy duration (DUR), and cumulative energy (ENER) for broadband synthetics between 0.1 and 10 Hz for the 2008 Mw5.4 Chino Hills earthquake. In this paper, our simulations are limited to low frequencies. We define the GOF score here as simply an equally weighted average of PGV, PGD, RS (1-10, 4-10, 2-4, and 1-2 s), FS (0.1-1, 0.1-0.25, 0.25-0.5, and 0.5-1.0 Hz), DUR, and ENER for the synthetics with respect to the observed seismograms filtered between 0.1 and 1.0 Hz, as already shown in Figs. 5, 6, and 7. Figure 10 shows the GOF scores calculated for the waveforms based on models K, AD, and AK at surrounding stations; the detailed scores are shown in Tables 1, 2, and 3 of the Appendix for each ground motion parameter to figure out the impacts in the synthetic seismograms. It is generally considered that the score is Bexcellent^for more than 80 %, Bgood^for more than 60 %, Bfair^for more than 40 %, and Bpoorf or the rest. As observed in the waveform comparison in the previous section, the GOF scores for both AD and AK models are generally better than those from the kinematic model (K). One of the reasons is that the simulations using model K were not computed using the 3D geological model. However, we think the principal reason is that the kinematic model may have smoothed the rupture process itself because of the way in which it was originally determined by inversion, e.g., constraints on the source time function and sub-fault size. The dynamic rupture model, on the other hand, can naturally describe any drastic change in the rupture process. We also note that at particular stations, the GOF scores are not good for any source model. For example, stations NIG010 and NIG011 are located relatively far from the source in the Niigata basin, where the surface waves become dominant and continue shaking until the end of the calculation (60 s). This is one reason why the score is not good. On the other hand, we know that the nearest station, NIG018 (in Kashiwazaki basin), is at a site that liquefied; it is not possible to compare the scores when the ground motion simulation is based on a linear-elastic calculation. However, these stations are also included in the averaged GOF scores, which are still in the Bfair^range for the dynamic source models AD and AK. Although the statistical analysis on the synthetics provides a means to quantify how good the results are among the different models, one should still try to understand the basic features of each seismogram. Fig. 10 Comparison of the GOF scores for the synthetics from the three source models (K, AD, and AK) with respect to the observed seismograms. The GOF score is an average of the six metrics (PGV, PGD, DUR, RS, FS, and ENER) calculated from the waveforms filtered between 0.1 and 1.0 Hz. Broken line shows the average GOF score of each source model over all the stations. Thick short line is the GOF score averaged for each station. GOF scales are from 0 (poor) to 100 (excellent). The detailed scores by each ground motion parameter are shown in Tables A1, A2 and A3 of the Appendix Kristekova et al. (2009) proposed a method to evaluate a GOF score in the time-frequency space in order to judge the local variation of the fit between the data and synthetic waveforms. Figure 11 shows the GOF score (normalized to 10) on the time-frequency envelope between the observed and synthetic ground velocities at station KSH-SG4 for all three components. Compared to the kinematic model K, the dynamic model AD for this station is slightly different for the EW component, but quite different for the NS component. This figure confirms that the dynamic model AK shows similar scores over the time-frequency range. The colored (yellow + orange) area of model AK is slightly smaller than in model K (by 23 % for the EW component). This comparison of models shows that the dynamic model can generate input ground motion equally well as the kinematic model. This suggests that a hybrid approach, i.e., combining a low-frequency dynamic rupture with any stochastic approach, could be used for generating the broadband ground motion in future applications, as suggested with various finite-source models (e.g., Graves and Pitarka 2015;Crempien and Archuleta 2015;Olsen and Takedatsu 2015).

Discussion and summary
The aim of this paper was to examine the performance of the dynamically simulated earthquake ruptures in computing ground motion for the purpose of quantitative seismic hazard analysis. For the 2007 Mw6.6 Chuetsu-oki, Japan, earthquake, we compute ground motion in a 3D geological model using a kinematic source model and two dynamic source models. To quantify the performance of the models, we apply a GOF criterion to the simulated ground motions. From the comparisons, we find that the kinematic model does not always have the best performance in reproducing the characteristics of the strong ground motion. In particular, it does not work well in the very near field and in the forward direction of the rupture propagation. This is probably because an abrupt change in the rupture Fig. 11 Time-frequency envelope GOF for the synthetics at SKH-SG4. The three source models are compared for three components (x: EW, y: NS, and z: UD) of the velocity seismograms. GOF scales are from 0 (poor) to 10 (excellent) process (rupture onset and changes in rupture velocity) may not be well simulated by the kinematic description. Two dynamic rupture models produce a similar ground motion radiation. However, the model with the asperities distributed on a single fault plane might be a projection of the rupture process on conjugate fault segments. The rupture process changes the wave radiation naturally due to the geometrical irregularities. From this perspective, the dynamic rupture model on a complex fault geometry produces a reasonable rupture scenario and wave radiation for practical applications.
For the 2007 Chuetsu-oki earthquake, the conjugate segments might be a reasonable causal source model, as inferred from the geodetic analyses (Nishimura et al. 2008;Aoki et al. 2008). Although the model parameters could be calibrated better, the rupture process on each segment could be very simple, represented by three phases: first, rupture on the NW-dipping segment; second, a dynamic rupture transfer between conjugate segments (which changes the wave radiation); and third, rupture on the SE-dipping segment. These features correspond to the three asperities commonly found from the seismological finite-source inversions. It is not possible to distinguish between the two rupture scenarios. To allow for possible mechanisms, one has to consider the possibility of rupture transfer from one segment to another and be included in probable rupture scenarios used to estimate the ground motion. Indeed, recent improvements of the geophysical observations often reveal complex fault geometries even for moderate-magnitude earthquakes, such as the 2009 Mw6.4 Suruga Bay, Japan (Aoi et al. 2010), and the Mw6.9 Iwate-Miyagi Nairiku, Japan, earthquakes (Fukuyama, 2015). In the seismological analyses, these earthquakes were mostly studied as ruptures on a single fault plane, an approximation good enough for regional or teleseismic scales.
There remains a scientific debate on the different aspects of the dynamic rupture process because the dynamics are technically difficult to solve and the frictional component of the faulting is difficult to study. However, progress over the last two decades helps us to understand various aspects of the dynamics of the earthquake mechanism. Dynamic rupture scenarios have been used to compute the ground motions for highseismic-hazard areas such as California (e.g., Olsen et al. 2009), Japan (Sekiguchi and Kase 2012), and Turkey (Aochi and Ulrich 2015). The parameter studies on dynamic rupture allow retrieving probable rupture scenarios and ground motions (e.g., Aochi and Ulrich 2015). Moreover, these dynamic simulations provide insight about the variability of the phenomena and extreme ground motions (e.g., Andrews et al. 2007).

Appendix
Tables 1, 2, and 3 summarize the GOF scores for each ground motion parameter by component for the three source models (K, AD, and AK), respectively. The average of the scores by each seismogram is shown in Fig. 10. DUR energy duration, PGV peak ground velocity, PGD peak ground displacement, RS response spectral acceleration (averaged for periods between 1 and 10 s (the entire period of study), 4 and 10, 2 and 4, and 1 and 2 s, respectively), FS Fourier spectrum (smoothed to reduce variance for frequencies between 0.1 and 1 Hz (the studied entire frequency), 0.1 and 0.25, 0.25 and 0.5, and 0.5 and 1 Hz, respectively), ENER cumulative energy (Olsen and Mayhew 2010). For some stations, we observe that the ground oscillation does not cease yet in Figs. 5, 6, and 7 by the end of the calculation (60 s). In such cases, we had to cut the end of DUR automatically by the end of the calculation. This may be one of the reasons why the DUR score is not high compared to the other parameters.