# On the effect of the 3-D regional geology on the seismic design of critical structures: the case of the Kashiwazaki-Kariwa Nuclear Power Plant

On the effect of the 3-D regional geology on the seismic design of critical structures: the case... Summary In this study, numerical investigation is performed on a realistic source-to-site earthquake scenario, with the aim to assess the role of complex 3-D geological structures on the predicted wavefield. With this respect, the paper pointedly targets the seismic response of nuclear power plants in near-field conditions and the verification of some simplified assumptions commonly adopted for earthquake ground motion prediction and site effects analysis. To this purpose, the Kashiwazaki-Kariwa Nuclear Power Plant (Japan) is assumed as reference case-study. In 2007, the nuclear site and its surroundings were struck by the Niigata-Ken Chūetsu-Oki seismic sequence, which caused some of the peak ground motion design limits to be largely overpassed. The dense observation network deployed at the site recorded a highly incoherent and impulsive earthquake ground motion. Many studies argued that the intricate syncline-anticline geology lying underneath the nuclear facility was highly responsible of the observed seismic response. Therefore, a physics-based numerical model of the epicentral area is built-up (≈60 km wide) and tested for small aftershocks, so to discount the effect of extended source on the synthetic site-response. The numerical model (based on the Spectral Element Method) reproduces the source-to-site wave propagation by embracing the effects of the surface topography along with the presence of the Japan Sea (i.e. the bathymetry, the coastline and the fluid–solid interaction). Broad-band (0–5 Hz) synthetic waveforms are obtained for two different aftershocks, located at the two opposite sides of the nuclear facility, aiming to assess the influence of the incidence angle the radiated wave field impinges the foldings beneath it. The effect of the folding presence is assessed by comparing it to a subhorizontally layered geology, in terms of numerical outcome, and by highlighting the differences with respect to the observations. The presence of an intricate geology effectively unveils the reason behind the observed ground motion spatial variability within a relatively small area, stressing its crucial role to properly reproduce the modification the wavefield undergoes during its propagation path towards the surface. The accuracy of the numerical exercise is discussed along with its results, to show the high-fidelity of these deterministic earthquake ground motion predictions. Computational seismology, Earthquake dynamics, Earthquake ground motions, Wave propagation 1 INTRODUCTION In the past two decades, the seismic hazard analysis and vulnerability assessment of large urban areas and critical structures and infrastructures have progressively taken advantage of the ever-increasing computational power available (Paolucci et al.2014). This outstanding technological and numerical progress seemingly broke through the evergreen and most stringent bottleneck in computational seismology: the impossibility to solve the complete source-to-site seismic wave propagation problem in a single-step analysis. This achievement paved the way to fully couple the large scale seismological models for the region of interest (i.e. the fault mechanism and the geological properties of the Earth’s crust), with local engineering models for geotechnical, site-effect and structural analyses, in the next few years. To this end, a plethora of extremely efficient numerical methods, for instance, the Finite Difference Method (FDM; Graves 1996; Pitarka et al.1998; Moczo et al.2007), the Finite Element Method (FEM; Taborda & Bielak 2011, 2013) and the Spectral Element Method (SEM; Faccioli et al.1997; Komatitsch & Vilotte 1998; Mazzieri et al.2011) has been progressively applied to the seismological problem. In spite of the inherent complexity and the huge dimensions of those computational models, their power is essentially embodied by the higher broad-band accuracy they provide (i.e. up to 4-5 Hz, De Martin 2011; Paolucci et al.2015), gradually bridging the gap between low-frequency source models obtained via waveform inversion techniques and the structural modal frequencies (i.e. up to 20 Hz). The mentioned enhancement in computational seismology was certainly granted by improved vectorization strategies to implement and tailor the mentioned numerical methods to massively parallel architectures (Dupros et al.2010; Mazzieri et al.2013; Göddeke et al.2014). Once shattered the computational barriers, a new holistic philosophy took place, driven by the deterministic modelling of the physics underlying each aspect of the earthquake phenomenon, for more accurate sensitivity analyses and uncertainty quantification of models and related parameters. All the ingredients (i.e. source, path and site effects) are naturally convolved in one-step all-embracing analysis, which is capable to predict a realistic seismic wavefield and to explain the observed time-histories in sedimentary deposits and/or at a global scale (De Martin 2011). A 20 yr history of 3-D physics-based simulations of past earthquakes has been written so far, with an overall increasing frequency upper limit (Graves 1996; Pitarka et al.1998; Bao et al.1998; Komatitsch et al.1999; Graves & Wald 2004; Day et al.2008; Gallovič et al.2010; Smerzini & Villani 2012; Aochi et al.2013; Paolucci et al.2015,among others). For a detailed excursus, please refer to Paolucci et al. (2015), who listed a selection of recent studies to produce earthquake ground shaking scenarios in large metropolitan areas. To clarify, by regional scale, it is hereafter assumed an average distance range between about 1 km (local scale) to 100 km (continental scale). Moreover, several numerical benchmarks have been recently performed during the past years, with the aim at testing the accuracy and the efficiency of different simulation methods when handling complex domain geometries and material interfaces. For instance, Chaljub et al. (2015) tested different 3-D numerical methods for earthquake ground shaking scenarios. They defined four canonical models inspired by the Mygdonian basin near Thessaloniki, Greece. The wider horizons opened by the described computational tools imply the need for detailed velocity structures, so to explain their impact on the radiated broad-band wavefield. Historically, long-period synthetic waveforms for large seismic events made use of the empirical Green’s function method, based on the self-similarity principle (Aki 1967) and on the representation theorem (Spudich & Archuleta 1987), and firstly proposed by Hartzell (1978) and Irikura (1983). Empirical approaches are however suitable just in the favourable yet rare condition of a vast number of recordings available for the site of interest or a very similar one (Chaljub et al.2015). However, a general concept was mediated from the mentioned traditional approach, that is, the idea that seismic observations recorded during small aftershocks identify with the system response to point-wise impulsive source excitation. From this point of view, the simulation of small aftershocks clarifies the impact of the propagating path and site effects, whose nature is manifold and sometimes counter-intuitive. One of the most investigated and debated aspects is the effect of the buried topography on the surface wavefield, due to the presence of complex geological interfaces such as basins, valleys, foldings and large heterogeneous soil deposits. The latter issues, along with the surface topography are responsible of wave scattering and diffraction, of body wave (de-)amplification, focusing and multiple reflection and sometimes of the presence of long-period surface waves. These effects may amplify due to proximity of the source to the site and, from a structural engineering point of view, they translate into a spatially incoherent incident motion, which entails the need for more refined seismic analysis when critical structures (such a nuclear power plants) are designed. For instance, Paolucci et al. (2015) recently provided an interesting insight on the role of the folded buried structure underneath the Po Plain, as one of the prominent causes of the ground motion distribution observed during the 2012 Emilia earthquake (Italy). Further interesting studies on the influence of the 3-D geological configurations on the characteristics of the ground shaking may be found in (Bard & Bouchon 1980; Faccioli & Vanini 2003; Gélis & Bonilla 2012; Chaljub et al.2015,among others). In this paper, the impact of a folded geology on the predicted ground shaking scenario of the site response of a nuclear power plant is assessed. The reference test case is represented by the 2007 MW6.6 Niigata-Ken Chūetsu-Oki earthquake (NCOEQ-2007), which struck a 100 km-wide area along the Japanese mid-west coastal area (depicted in Fig. 1), encompassing the Kashiwazaki-Kariwa Nuclear Power Plant (KKNPP). Figure 1. View largeDownload slide Map of the Niigata region, surrounding the Kashiwazaki-Kariwa Nuclear Power Plant (black square). The NCOEQ-2007 main shock is indicated as a red star, whereas coloured squares indicate the relocated aftershocks positions (TEPCO 2007). ASP1, ASP2 and ASP3 indicate the three major asperities deducted by Shiba (2008), by means of waveform inversion. Figure 1. View largeDownload slide Map of the Niigata region, surrounding the Kashiwazaki-Kariwa Nuclear Power Plant (black square). The NCOEQ-2007 main shock is indicated as a red star, whereas coloured squares indicate the relocated aftershocks positions (TEPCO 2007). ASP1, ASP2 and ASP3 indicate the three major asperities deducted by Shiba (2008), by means of waveform inversion. The KKNPP is the largest nuclear power plant in the world, with its seven units deployed around a 1.5 km wide area. The facility is located in the southern Niigata prefecture, adjacently the affected area (the Joyner-Boore and rupture distances were RJB = 0 km and RRUP = 16 km respectively; Yee et al.2011) and specifically on the hanging wall of the buried reverse fault which caused the earthquake. It has been reported that the site experienced nearly twice the value of Peak Ground Acceleration (PGA) assumed for the plant design (Pavlenko & Irikura 2012), but even despite the large ground settlement occurred, the redundant and over-designed structures incurred into minor damage (Kayen et al.2007). Nevertheless, in this paper, the spatial variability of the ground motion observed within the site area is investigated, based on previous studies which attributed it to the presence of a geological folding structure near the KKNPP (Tokumitsu et al.2009; Watanabe et al.2009; Hijikata et al.2010). For instance, Watanabe et al. (2009) and Watanabe et al. (2011) performed some 2-D FEM analyses, reliable up to 5.0 Hz, considering two cross-sectional of the mentioned folded structure. They investigated and successfully reproduced the amplification factors at two sites within the KKNPP, picturing the higher amplification occurred in the southernmost area. Following their work, Tsuda et al. (2011) extruded the 2-D cross sections along its directrix (approximately following the NW-SE coastline) and performed an FDM simulation of the NCOEQ-2007 main shock and three aftershocks, depicting the regional wavefield and the effect of an extended fault rupture, up to a frequency of 4.0 Hz. Aochi et al. (2013) tested three different 3-D geological structures (by simulating two aftershocks in the NCOEQ-2007 cluster) in an FDM simulation to reproduce the main shock scenario at long periods (i.e. up to 0.5 Hz). Aochi & Yoshimi (2016) improved the latter computational model, including the best available geological model and testing two different source models (both kinematic and dynamic descriptions, in the 0.1–1.0 Hz frequency band). Further pieces of research concerning the 3-D numerical simulations of the NCOEQ-2007 earthquake are the one provided by Kawabe & Kamae (2010) and Quinay et al. (2013). The formers simulated the wave propagation in a 3-D structure model, in a frequency range 0.05–1.6 Hz. Quinay et al. (2013) performed instead an entire fault-to-structure analysis, including the crustal geology, the soil strata and the KKNPP buildings. They tackled and solved this cumbersome computational problem by means of a two-step multiscale analysis based on (1) a solution of a low-resolution model (at a regional scale), used as an input boundary condition to the higher-resolution one (at an engineering-length scale). The computational scheme (formulated by means of the FEM) implied a domain decomposition with pre-partitioning, solved on parallel architecture. The authors obtained realistic results, compared to the observations in the Niigata region and at the KKNPP, being restrained to the frequency domain 0.1–1.0 Hz. In this context, most of the studies targeting the regional wavefield and the seismic source characterization lack of poor frequency resolution (see for instance Aochi & Yoshimi 2016; Quinay et al.2013), whereas models with higher frequency ranges (such as Watanabe et al.2009; Tsuda et al.2011) are focalized on the KKNPP subsurface structure. In this work, a 3-D regional-scale earthquake scenario of the area was built-up for two of the aftershocks belonging to the NCOEQ-2007 sequence and indicated as blue and green stars respectively in Fig. 1. Section 2 describes the main aspects of the observed waveforms at the KKNPP, focusing on the impulsive nature of the recorded velocigrams and the incoherent amplification between the southwest and northeast part of the site, placed at approximately 1 km from each other. The details on the geological structure in the Niigata region are provided in Section 3, stressing the methodological approach followed to conceive a fancy geological model suitable for reliable broad-band earthquake ground-motion simulations. The selection of the coarser and horizontally layered geology is described in Section 4 (the so called model LARGE), with the aim at reproducing the regional wavefield (solved by means of the Spectral Element Method, SEM) in a rather low-frequency range (up to fmax = 3.75 Hz, given the lack of a more detailed information at this scale) but constrained by several recording stations in the Niigata region. The impact of the folding is assessed in a second phase (Section 5) by plugging its 3-D configuration into the previous regional model (and obtaining the SEM model called SMALL). The outcome of numerical modelling from two aftershocks occurred at opposite sides of the KKNPP is analysed, comparing the synthetics (filtered at 5.0 Hz) at different locations within the nuclear plant. Special attention is given to the differences in terms of response spectra (Sa) at the ground surface, highlighting the discrepancies at different locations (corresponding to different Units of the nuclear site) and to their variation according to the aftershock position with respect to KKNPP. The comparison with the observations stress the need for a more accurate prediction of incident wave motion to be used as input for structural analyses, dropping the simplified assumption of layered half-space and supporting the importance of the complex geology as responsible of the ground motion incoherence at surface. 2 THE NCO EARTHQUAKE: OBSERVED STRONG GROUND MOTION CHARACTERISTICS The KKNPP site consists of seven Units, grouped into two blocks: (1) Units 1–4, built in the southwest part of the site and (2) Units 5–7 placed approximately 1.5 km away, in the northeast corner. Fig. 2 shows a sketch of the plant, highlighting the recording devices downhole and at surface. Devices 1G1 (green) and 5G1 (orange) both belong to the most recent (2004) recording network, and they are placed at the surface (G.L. 0 m), respectively close by Unit 1 and Unit 5–7. On the other hand, four devices (SG1–SG4) are deployed along the vertical array KSH (blue) located at the KKNPP Service Hall, up to a depth of 250 m; other five devices are deployed downhole close by Unit 5 (KK5, red), till 312 m of depth; four devices are finally located in the surroundings of Unit 1 (KK1, magenta) with devices G07–G10 reaching 255 m of depth. Figure 2. View largeDownload slide Map of the KKNPP site (courtesy of TEPCO 2007). The five coloured squares indicate the three recording stations downhole (KSH at Service Hall (blue), KK1 at Unit 1 (magenta), KK5 at Unit 5 (red)) and the two surface ones free-field (1G1 at Unit 1 (green) and 5G1 at Unit 5 (orange)). The KSH station downhole, along with recording devices at surface 1G1 and 5G1 entirely recorded the NCOEQ-2007 main shock. The devices were oriented with respect to plant North (NS-KKNPP), which differs from the real geodetic North of an angle θ = 18° 54$$^{^{\prime }}$$ 51$$^{^{\prime \prime }}$$. Moreover, some of them suffered of an azimuthal deviation, whose value was provided by TEPCO (2007). Figure 2. View largeDownload slide Map of the KKNPP site (courtesy of TEPCO 2007). The five coloured squares indicate the three recording stations downhole (KSH at Service Hall (blue), KK1 at Unit 1 (magenta), KK5 at Unit 5 (red)) and the two surface ones free-field (1G1 at Unit 1 (green) and 5G1 at Unit 5 (orange)). The KSH station downhole, along with recording devices at surface 1G1 and 5G1 entirely recorded the NCOEQ-2007 main shock. The devices were oriented with respect to plant North (NS-KKNPP), which differs from the real geodetic North of an angle θ = 18° 54$$^{^{\prime }}$$ 51$$^{^{\prime \prime }}$$. Moreover, some of them suffered of an azimuthal deviation, whose value was provided by TEPCO (2007). The three vertical arrays KSH, KK5 and KK1 belong to the oldest recording network, where surface devices 1G1 (Unit 1, G.L. 0 m) and 5G1 (Unit 5, G.L. 0 m) are deployed at surface. It is worth noting that just one (KSH, Service Hall) out of three vertical array entirely recorded the main shock, along with 1G1 and 5G1. According to Watanabe et al. (2009), the strong motion records at Unit 1 were significantly larger than those in the surrounding of Unit 5. Many authors (e.g. Uetake et al.2008; Miyake et al.2010; Tsuda et al.2011) stressed the impulsive nature of the recorded time-histories. Large velocity pulses are usually observed nearby the fault rupture and coupled with their corresponding large peak displacement, they considerably enhance the structural damage potential (Cox & Ashford 2002). In the NCOEQ-2007 case, Uetake et al. (2008) first noticed the three significant pulses at KSH array, associating them to the three major asperities identified (via waveform inversion) on the fault plane. As a matter of fact, velocity pulses may be effectively extracted from the time-histories at KSH (SG1, G.L. −2.4 m and SG4, G.L. −250 m), 1G1 and 5G1 (i.e. the only devices which entirely recorded the main shock), by employing the ranking criterion proposed by Baker (2007) (excluding late arrivals and small events). Fig. 3 shows an example of the mentioned pulse-like nature of the velocigrams recorded. Figure 3. View largeDownload slide Comparison between velocity signal (black traces) and extracted pulses (red traces), according to Baker (2007), for devices (a) SG1 (G.L. −2.4 m), (b) SG4 (G.L. −250 m) in the KSH vertical array (Service Hall) and for (c) 1G1 (G.L. 0 m) and (d) 5G1 (G.L. 0 m). KKZ indicates the set of new seismometers, installed in 2004, to which 1G1 and 5G1 belong. The in plane azimuthal angle θ at which those pulses were extracted is indicated above the velocigrams. Figure 3. View largeDownload slide Comparison between velocity signal (black traces) and extracted pulses (red traces), according to Baker (2007), for devices (a) SG1 (G.L. −2.4 m), (b) SG4 (G.L. −250 m) in the KSH vertical array (Service Hall) and for (c) 1G1 (G.L. 0 m) and (d) 5G1 (G.L. 0 m). KKZ indicates the set of new seismometers, installed in 2004, to which 1G1 and 5G1 belong. The in plane azimuthal angle θ at which those pulses were extracted is indicated above the velocigrams. The average pulse period ranges around 2.5 s for SG1, SG4 and 5G1, whereas 1G1 identifies with a shorter period pulse (1.5 s). Remarkable differences (in terms of amplitude, period and phase) are evident at a first glance, seemingly indicating an incoherency of the ground motion all around the nuclear facility zone. The observed site response at the surface is definitely more intense at the Service Hall and at Unit 1. It is worth noting that the velocity pulses were identified and extracted at different azimuthal directions (with respect to the real North) for SG1 (KSH, surface) SG4 (KSH, depth), 1G1 (Unit 1, surface) and 5G1 (Unit 5, surface) respectively, the polarization of the impulsive incident wavefield spatially varies at the site scale, along its travelling path towards the surface. Another interesting aspect resides in the fact that later pulses, whose nature has been correlated to the progressive fault rupture of the three asperities, may be identified at each location (Uetake et al.2008; Miyake et al.2010). Tsuda et al. (2011) observed that the first two pulses were rather comparable in terms of maximum amplitude, whereas the third one (allegedly coming from the third asperity ASP3, represented as green rectangle in Fig. 1) is different for two sites: the peak ground velocity on the south side of the KKNPP (KK1, downhole array at Unit 1) is much larger than that on the north side (KK5, downhole array at Unit 5) (Tsuda et al.2011). 3 CONSTRUCTION OF THE HYBRID GEOLOGICAL MODEL The subsurface geological structure underlying the Niigata region has been proven to be rather intricate. Several 1-D and 3-D models of the velocity structure have been proposed in the past, based on geological and geophysical explorations (see, for instance, Shinohara et al.2008). The latter set of geological configurations have been tested and tuned by several authors, by exploiting the long aftershock sequence occurred in the main shock aftermath and then exploited to simulate the main shock itself. Table 1 resumes the ensemble of ground shaking analyses performed on the KKNPP site, during the NCOEQ-2007. Kawabe & Kamae (2010) simulated the wave propagation in a 3-D structure model provided by the Japan Nuclear Energy Safety Organization (JNES, internal report, 2005) in a frequency range 0.05–1.6 Hz. Aochi et al. (2013) and Aochi & Yoshimi (2016) performed a complete FDM analysis of the NCOEQ-2007 scenario, comparing three regional velocity models: (1) the model proposed by the Earthquake Research Institute (ERI; University of Tokyo) obtained by P- and S-wave traveltime double-difference tomography on the data provided by Shinohara et al. (2008); Kato et al. (2008, 2009), with a grid resolution of 3 km × 5 km × 3 km (vertical) and with a minimum shear wave velocity VS of 866 m s−1; (2) the model proposed by the National Research Institute for Earth Science and Disaster Resilience (NIED), taken from Fujiwara et al. (2009) from the seismic reflection results and observed H/V spectra (available on the Japan Seismic Hazard Information Station (J-SHIS), website: http://www.j-shis.bosai.go.jp), with a 1 km × 1 km resolution and minimum VS equal to 350 m s−1; (3) the model proposed by the Geophysical Survey of Japan (GSJ) and found in Sekiguchi et al. (2009), which is an improvement of the NIED velocity structure for the Niigata area (the resolution is 0.5 km × 0.5 km), carefully tuned upon the regional variation of material parameters (VSmin = 400 m s−1). The authors concluded that the ERI model was less preferable than the other two in order to reproduce the regional wave field, mainly due to the poor resolution of the shallow deposits and despite the low frequency range of the analyses. Moreover, they argue the use of the JNES model (Kawabe & Kamae 2010) for the regional-wise purposes, being however more precise around KKNPP. The crust material, topography, and subsurface data provided by the J-SHIS were used by Quinay et al. (2013) to carry out a complete 3-D analysis up to 1.0 Hz, studying the effect of the synthetic wave field on the structural behaviour of one of KKNPP reactor building. Aochi & Yoshimi (2016) improved the previous results by Aochi et al. (2013) by enlarging the frequency band up to 1.0 Hz and employing the GSJ model, although their study targeted mostly the effect of different seismic source models and numerical implementation. Table 1. Summary of ground shaking earthquake ground shaking scenarios of the NCOEQ-2007, targeting the KKNPP. FDM, Finite Difference Method; FEM, Finite Element Method; SEM, Spectral Element Method. Fold, Tokumitsu et al. (2009); JNES-2005, internal report; ERI, Shinohara et al. (2008); Kato et al. (2008, 2009); JNES-2008, JNES (2008); NIED, Fujiwara et al. (2009); GSJ, Sekiguchi et al. (2009); Reference  Dimension  Geology  Method  Model Size [km]  fmax [Hz]  min(VS) $$[\rm {{m}\,{s}}^{-1}]$$  Watanabe et al. (2009)  2-D  Fold  FEM  7.6 × 4.8  5.0  700  Kawabe & Kamae (2010)  3-D  JNES-2005  ?  ?  1.6  700  Ducellier & Aochi (2010)  3-D  ERI  FDM  110 × 120 × 30  0.866  866  Tsuda et al. (2011)  3-D  JNES-2008 + Fold  FDM  50 × 50 × 20  4.0  700  Aochi et al. (2013)  3-D  ERI/NIED/GSJ  FDM  110 × 120 × 30  0.5  866/350/400  Quinay et al. (2013)  3-D  NIED  FEM  110 × 120 × 30  1.0  866/350/400  Aochi & Yoshimi (2016)  3-D  GSJ  FDM  110 × 120 × 30  1.0  400  LARGE (this work)  3-D  Aochi2013  SEM  90 × 83 × 82  3.75  1020  SMALL (this work)  3-D  Aochi2013+Fold  SEM  68 × 50.8 × 50.8  5.0  700  Reference  Dimension  Geology  Method  Model Size [km]  fmax [Hz]  min(VS) $$[\rm {{m}\,{s}}^{-1}]$$  Watanabe et al. (2009)  2-D  Fold  FEM  7.6 × 4.8  5.0  700  Kawabe & Kamae (2010)  3-D  JNES-2005  ?  ?  1.6  700  Ducellier & Aochi (2010)  3-D  ERI  FDM  110 × 120 × 30  0.866  866  Tsuda et al. (2011)  3-D  JNES-2008 + Fold  FDM  50 × 50 × 20  4.0  700  Aochi et al. (2013)  3-D  ERI/NIED/GSJ  FDM  110 × 120 × 30  0.5  866/350/400  Quinay et al. (2013)  3-D  NIED  FEM  110 × 120 × 30  1.0  866/350/400  Aochi & Yoshimi (2016)  3-D  GSJ  FDM  110 × 120 × 30  1.0  400  LARGE (this work)  3-D  Aochi2013  SEM  90 × 83 × 82  3.75  1020  SMALL (this work)  3-D  Aochi2013+Fold  SEM  68 × 50.8 × 50.8  5.0  700  View Large 3.1 Description of the layered geology of the Niigata region As a general remark, the lowest shear-wave velocity varies from 350 m s−1 to 866 m s−1, for the studies previously, estimated for low-frequency analyses (i.e. ground motion reproduction and a few waveform inversions of the source mechanism), and eventually based on the interpretation of several 1-D/2-D models, which Aochi et al. (2013) proved to perform fairly well in reproducing the complex wavefield radiated, although phase shifts must be taken into account, due to modified traveltimes. Fig. 4 shows some of the mentioned layered velocity structures available for the Niigata region. D&A2010-1 and D&A2010-2 were proposed by Ducellier & Aochi (2010), in their numerical simulations performed for the NCOEQ-2007 main shock. With the same purposes, Aochi et al. (2013) provided a simplified 1-D profile, tagged as Aochi2013. Cire.2008 velocity values were employed by Cirella et al. (2008) to perform a waveform inversion of the NCOEQ-2007 source mechanism, whereas Shino.2008 was obtained by exploiting the data of the seismic survey conducted by Shinohara et al. (2008) and used to relocate the aftershock sequence. Finally, profile indicated as NIG004, NIG016, NIG026 and NIGH12 refer to the namesake K-NET and KiK-net stations. Figure 4. View largeDownload slide Set of regional 1-D velocity structures (VP and VS) for the Niigata region, estimated by several authors. D&A2010-1 and D&A2010-2 were proposed by Ducellier & Aochi (2010). Aochi et al. (2013) and Cirella et al. (2008) provided the profiles tagged as Aochi2013 and Cire.2008 respectively. Shino.2008 was obtained by Shinohara et al. (2008) and profiles indicated as NIG004, NIG016, NIG026 and NIGH12 refer to the namesake K-NET and KiK-net stations (the NIED strong ground motion networks, covering the Japan territory) and they were estimated by K. Koketsu (source: http://taro.eri.u-tokyo.ac.jp/saigai/chuetsuoki/source/index.html) Figure 4. View largeDownload slide Set of regional 1-D velocity structures (VP and VS) for the Niigata region, estimated by several authors. D&A2010-1 and D&A2010-2 were proposed by Ducellier & Aochi (2010). Aochi et al. (2013) and Cirella et al. (2008) provided the profiles tagged as Aochi2013 and Cire.2008 respectively. Shino.2008 was obtained by Shinohara et al. (2008) and profiles indicated as NIG004, NIG016, NIG026 and NIGH12 refer to the namesake K-NET and KiK-net stations (the NIED strong ground motion networks, covering the Japan territory) and they were estimated by K. Koketsu (source: http://taro.eri.u-tokyo.ac.jp/saigai/chuetsuoki/source/index.html) 3.2 Description of the local 3-D geology Despite the good approximation of the regional incident wave motion obtained with simplified 1-D and 3-D geological configurations in the surroundings of the epicentral area, there is consensus among many researchers (Uetake et al.2008; Tokumitsu et al.2009; Watanabe et al.2009; Hijikata et al.2010) that none of them is suitable to picture the spatial variability of the ground motion observed at KKNPP. With this regard, Uetake et al. (2008) and Tokumitsu et al. (2009) related this aspect to the presence of a folded structure underneath KKNPP (characterized via boring and seismic reflection surveys, Kobayashi et al.1995) and composed by the Madonosaka syncline, interposed between the Ushirodani and Chuo-Yutai anticlines (Gürpinar et al.2017). The folding cross-section extends for approximately 7 km across the Japanese coastline, up to 2.5 km down depth and with its hinge axes (i.e. the directrixes) strike at N55°E, as shown in Figs 5(a) and (b). The geographical distribution of this peculiar sediment conformation shears through the nuclear facility, with the Ushirodani anticline and Madonosaka syncline placed below Unit 5 and Unit 1 respectively. The Chuo-Yatai anticline is indicated as well, which seemingly passes in the Service Hall surroundings (array KSH). A detailed description of the seven strata defining the folding is provided in Section 5. Watanabe et al. (2009) constructed a 2-D FEM model of the syncline-anticline conformation, from boring and seismic reflection survey. Simulation results show good agreement with the observed strong motion records at Unit 1. Tsuda et al. (2011) constructed a 3-D model by taking into consideration seven 2-D sections across the folding area, and interpolating them. The joined up the folding model with the regional one (the 3-D velocity model built by the Japan Nuclear Energy Safety Organization (JNES), which covers broad area including Chuetsu area (JNES 2008). 3.3 Proposed hybrid model of the KKNPP earthquake scenario Given the complex geological structure described so far, and with the twofold aim to reproduce the overall regional wavefield, yet focusing on a realistic broad-band simulation of the incident ground motion in the KKNPP surroundings, a hybrid geological model of the Niigata region was herein constructed in two steps, progressively adding all the available information. First of all, the layered Aochi2013 profile (described in Section 3.1) was selected as representative for the propagation of the long period ground motion (i.e. 0.1–3.75 Hz) at a regional scale (i.e. approximately 90 km × 83 km × 82 km). From now on, it will be referred as to the LARGE model. Its considerable extension was mainly chosen to grant a minimal coverage of the NW regional wavefield. Thus LARGE model was designed so as to include part of the Sado Island in front of the Niigata region, where the K-NET station NIG004 is located. The Japan sea presence was disregarded at this step, along with related bathymetry and coastline (see Fig. 5b). Figure 5. View largeDownload slide (a) Schematic structural map of the series of Ushirodani anticline—Madonosaka syncline, located underneath KKNPP. Red, orange and fuchsia points indicate the KKNPP site and AS1 and AS2 epicentres, respectively. The white line represents section SC. (b) Sketch of the extension of the computational models (LARGE and SMALL) considered in this study. Figure 5. View largeDownload slide (a) Schematic structural map of the series of Ushirodani anticline—Madonosaka syncline, located underneath KKNPP. Red, orange and fuchsia points indicate the KKNPP site and AS1 and AS2 epicentres, respectively. The white line represents section SC. (b) Sketch of the extension of the computational models (LARGE and SMALL) considered in this study. Section 4 presents the calibration of LARGE model, dealing with the choice of an appropriate source time function and topographical effect. The LARGE model was trimmed and adjusted so as to plug the 3-D Ushirodani anticline—Madonosaka syncline—Chuo-Yatai anticline structure (widely described in Section 3.2) into it, following the previous work of Tsuda et al. (2011) (see Section 5). The 1-D profile, issued from Aochi2013, was welded to the syncline-anticline system by adjusting and linearly smoothing the uppermost crustal material discontinuities to reconnect with the external layered geology of LARGE model. 4 SIMULATION OF REGIONAL WAVEFIELD: LARGE MODEL 4.1 Preliminary selection of a suitable geological profile The epistemic uncertainty on the regional geology was outlined and solved by first employing the Wave-Number-Integration method (WNI), proposed and implemented by Hisada (1994, 1995, 2008), that is, a semi-analytical method that simulates the complete 3-D wave propagation field (based on the computation of static and dynamic Green’s functions of displacements and stresses) radiated from an extended kinematic seismic source in a viscoelastic horizontally layered half-space. The MJMA4.4 aftershock of July 16, 21.08h (JST), was considered, which nucleated near-by one of the three major fault asperities at a depth of 11 km (Tsuda et al.2011) determined by detecting the polarity for the P-wave first motion (performed by the Hi-net (Obara et al.2005)) and by considering the Centroid Moment Tensor analysis provided by F-net data (Kubo et al.2002). The KiK-net and K-NET recording stations were employed at this point, along with the recorded time-histories at the KKNPP (device KSH-SG4, placed at G.L.-250m at the Service Hall of the nuclear site, see Fig. 2). Several WNI analyses were performed, testing the estimated geology in Fig. 4. This preliminary investigation ended up with the choice of the Aochi2013 soil profile proposed by Aochi et al. (2013) and whose velocity and attenuation properties are summarized in Table 2. Fig. 6 shows a fairly good comparison obtained at several stations in the surroundings of the KKNPP. In accordance with Aochi et al. (2013), the recordings and the synthetics were band-pass filtered between 0.1–0.5 Hz and aligned at the time-step corresponding to the 1 per cent of the respective Arias intensity. Although a satisfactory match has been obtained for other soil profiles (in fact, similar results were obtained for the crustal models tagged as Cire.2008 and Shin.2008 in Fig. 4), Aochi et al. (2013) introduced a sequence of softer layers (Fig. 4) at shallow depths, with respect to the crustal models estimated for the K-NET/KiK-net stations (NIG004, NIG016, NIG026, NIGH12). Despite the local alignment, some phase shift among the simulations and the recordings is still visible for NIG004 and NIGH12, whereas in the near field (cf. NIG018, NIG016) the fit is quite good (probably due to the simple propagation path, Aochi et al.2013)). Finally, from Fig. 6, it appears that the simulations performed at KKNPP site are satisfactory only along the NS direction, whereas synthetics are de-amplified along the EW direction. This might be an effect of the shallow local geology of the Niigata basin. With that being said, it is well known that the choice of a 1-D velocity structure tends to accentuate the ground motion directionality. Figure 6. View largeDownload slide Simulations performed by WNI (red line) method of the NCOEQ-2007 MJMA4.4 aftershock of July 16, 21.08h (JST). Blue waveforms (REC) represent the recorded velocigrams, red waveforms (WNI) the synthetic ones. Both records and synthetics were base-line corrected and bandpass filtered between 0.1 and 0.5 Hz. Synthetic waveforms were obtained by considering the soil profile Aochi2013 (Fig. 4). Velocigrams (in cm s−1) are herein magnified by a factor 1000 and aligned with respect to the 1 per cent of their Arias intensity. Figure 6. View largeDownload slide Simulations performed by WNI (red line) method of the NCOEQ-2007 MJMA4.4 aftershock of July 16, 21.08h (JST). Blue waveforms (REC) represent the recorded velocigrams, red waveforms (WNI) the synthetic ones. Both records and synthetics were base-line corrected and bandpass filtered between 0.1 and 0.5 Hz. Synthetic waveforms were obtained by considering the soil profile Aochi2013 (Fig. 4). Velocigrams (in cm s−1) are herein magnified by a factor 1000 and aligned with respect to the 1 per cent of their Arias intensity. Table 2. Properties of the Aochi2013 profile (Aochi et al.2013). z represents the depth of the upper layer surface, VP and VS are the P-wave and S-wave velocities respectively, QP and QS the quality factors for P-wave and S-wave respectively. The ⋆ indicates the interface chosen to plug the folding structure into the 1-D Aochi2013 profile, granting a smooth transition from one model to the other. z[km]  $$V_P [\rm {{km}\,{s}}^{-1}]$$  QP[1]  $${V}_{S} [\rm {{km}\,{s}}^{-1}]$$  QS[1]  0.0  2.28  200.0  1.02  100.0  0.5  2.57  228.7  1.23  100.0  1.0  2.93  257.1  1.51  100.0  1.5  3.09  293.3  1.63  100.0  2.0  3.25  309.1  2.09  100.0  3.0  3.68  325.9  2.33  100.0  4.0  4.03  368.2  2.49  100.0  ⋆5.0  4.30  403.2  2.75  100.0  6.0  4.55  430.9  2.89  100.0  7.0  4.76  455.8  3.10  100.0  8.0  5.00  476.4  3.40  100.0  9.0  5.37  500.2  3.76  100.0  10.0  5.88  537.5  3.80  100.0  14.0  6.51  588.7  3.85  250.0  z[km]  $$V_P [\rm {{km}\,{s}}^{-1}]$$  QP[1]  $${V}_{S} [\rm {{km}\,{s}}^{-1}]$$  QS[1]  0.0  2.28  200.0  1.02  100.0  0.5  2.57  228.7  1.23  100.0  1.0  2.93  257.1  1.51  100.0  1.5  3.09  293.3  1.63  100.0  2.0  3.25  309.1  2.09  100.0  3.0  3.68  325.9  2.33  100.0  4.0  4.03  368.2  2.49  100.0  ⋆5.0  4.30  403.2  2.75  100.0  6.0  4.55  430.9  2.89  100.0  7.0  4.76  455.8  3.10  100.0  8.0  5.00  476.4  3.40  100.0  9.0  5.37  500.2  3.76  100.0  10.0  5.88  537.5  3.80  100.0  14.0  6.51  588.7  3.85  250.0  View Large 4.2 Verification of the model upgrade (WNI→SEM) A Spectral Element Method software, SEM3D, was used to solve the 3-D wave propagation problem so as to verify the code with the WNI semi-analytical model tuned in the earlier stages of the parametric analysis. The chosen 1-D geological model was implemented in 3-D large scale SEM numerical model. No topography was introduced at first, so as to compare the synthetics with WNI results. Five Gauss–Lobatto–Legendre (GLL) integration nodes (the integration points for the SEM formulation) were chosen to grant the minimum wavelength (≈500 m) a satisfactory discretization (i.e. to reach a frequency upper bound of ≈2 Hz, with a third order polynomial degree). In the SEM formulation, material properties are assigned to each GLL point. Therefore, two strategies were investigated to feature the numerical model with the varying material properties, namely (1) a direct meshing approach that adapts the computational grid to the layer-to-layer interface, so to coherently reproduce the abrupt impedance contrast along it and (2) a not honouring approach that associate to the whole independently meshed domain a space-varying heterogeneous material map (in this case, transverse isotropic in the horizontal plane). In this standard test case (i.e. a flat and bounded horizontally layered half-space) the second strategy does not seem particularly appealing since the material interfaces are planar and horizontal and a structured mesh was adopted. The numerical verification of the model upgrade (WNI→SEM3D) was satisfactory. For instance, Fig. 7 shows a synoptic comparison between synthetic waveforms filtered between 0.1 and 1.0 Hz. Although no computational gain is achieved by using one or the other approach at this early stage, the verification test was performed in view of implementing the not honouring approach in the non-structured mesh domain, extruded at depth from the local Digital Elevation Model. Figure 7. View largeDownload slide Numerical verification of SEM3D code, against the semi-analytical solution obtained via WNI method (0.1-1.0 Hz). Synthetic velocity waveforms obtained at (a) the KKNPP site (KSH-SG4, G.L. −250 m, focal distance 15.1 km) and at (b) the K-NET station NIG018 (focal distance of 20.1 km) are portrayed. A good match is achieved between the WNI method (WNI, red line) and SEM3D numerical model equipped either with a layered geology (LAY, blue line) or with a not honouring approach (GRD, green line) Figure 7. View largeDownload slide Numerical verification of SEM3D code, against the semi-analytical solution obtained via WNI method (0.1-1.0 Hz). Synthetic velocity waveforms obtained at (a) the KKNPP site (KSH-SG4, G.L. −250 m, focal distance 15.1 km) and at (b) the K-NET station NIG018 (focal distance of 20.1 km) are portrayed. A good match is achieved between the WNI method (WNI, red line) and SEM3D numerical model equipped either with a layered geology (LAY, blue line) or with a not honouring approach (GRD, green line) 4.3 Effect of the surface topography A further step of the numerical exercise consisted into adding the surface topography (TOPO) to the previously calibrated flat model (FLAT). 15 Gauss-Lobatto-Legendre (GLL) integration nodes were chosen to propagate a seismic wavefield up to ≈3.75 Hz (considering, as reference, a third-order polynomial interpolation, Paolucci et al.2015). The shear-wave velocity profile (Aochi2013 model) was implemented via a not honouring approach. The Sea of Japan is not taken into account: a flat solid surface was considered instead. Fig. 8 shows the comparison between the accelerograms obtained at the ground level at KKNPP site, for the FLAT (red line) and TOPO (blue line) models. The waveforms were band-pass filtered between 0.05 and 3.75 Hz. Early wave arrivals at higher frequencies are observed when the surface topography is included in the model, although it does not seem to consistently affect the radiated wavefield. This aspect might be due to the fact that the Japan sea was disregarded. Figure 8. View largeDownload slide Numerical verification of SEM3D code: effect of the surface topography. Synthetic waveforms obtained at KKNPP (at surface) are portrayed: (a) EW component, (b) NS component. Both synthetics were compared in the frequency band 0.05–3.75 Hz. Figure 8. View largeDownload slide Numerical verification of SEM3D code: effect of the surface topography. Synthetic waveforms obtained at KKNPP (at surface) are portrayed: (a) EW component, (b) NS component. Both synthetics were compared in the frequency band 0.05–3.75 Hz. 4.4 Comparison with ground motion prediction equations The a priori selection of a suitable Source Time Function (STF) for a point-wise double-couple approximation of the seismic source is not obvious (Bizzarri 2014). To clarify this fundamental aspect, a brief parametric analysis addressing different STFs and the respective featuring parameters was performed. For instance, Aochi et al. (2013) simulated two aftershocks as double-couple point-wise sources, with smooth bell-shaped STF (i.e. cubic B-spline function) of a 0.5 s duration at the hypocentre. Specifically, two classical STF are addressed, namely (1) the so called Bouchon’s ramp (or smoothed ramp function) $$u_S^{\text{TNH}}\left(t\right)$$ and (2) the source model proposed by Brune (1970) $$u_S^{\text{EXP}}\left(t\right)$$. The former is well suited to approximate long-tailed Moment Rate Functions (MRF; Duputel et al.2013). Two are the free parameters that tune the shape and the position in time of the two mentioned STFs: the maximum displacement offset umax and the generic rise-time τR. The maximum slip value umax is spontaneously determined by relating it to the scalar seismic moment value (for the MJMA4.4 aftershock event, the F-net indicated a seismic scalar moment M0 = 5.21× 1015 Nm), whereas τR is rather delicate to be estimated. Although no explicit definition of the rise-time is available neither for $$u_S^{\text{TNH}}\left(t\right)$$ nor for $$u_S^{\text{EXP}}\left(t\right)$$ (Bizzarri 2014), τR represents the duration of the slip rate pulse (Bizzarri 2014; De Martin et al.2007). Whilst estimated fairly easily for strong ground shaking, empirical correlations may provide longer rise-times, due to the fact that they are based upon registrations of larger earthquake events, polluted by the effect of the rupture time along the fault plane. To assess the influence of τR on the radiated regional-wave field, three sources were tested on the LARGE model (reported in Table 3). Two values of τR were chosen, being τR1 the value of the rise-time proposed in Dreger et al. (2007), estimated as $$\tau _{R1} = 10^{0.5\left({\rm M}_{\rm W}-6.69\right)}$$, and τR2 being derived from the expression of the semi-duration proposed by Duputel et al. (2012) $$\tau _{R2}/2 = 1.2\cdot 10^{-8} \root 3 \of {{\rm M}_0\left[{\rm dyne.cm}\right]}$$. The latter value is closer to the rise-time indicated by Tsuda et al. (2011), based on the width of S-wave pulse of the observed waveforms at rock sites of F-net stations adjacent to the source (and equal to 0.7 s). Table 3. Source Time Functions (STFs) chosen to represent the point-wise seismic moment time-history. Brune’s pulse $$u_S^{{\rm EXP}}$$ refers to the STF proposed by Brune (1970), whereas the Bouchon’s ramp $$u_S^{{\rm TNH}}$$ refers to the STF proposed by Bouchon (1981). H(t) represent the Heaviside function. STF  Function  $$u_S^{{\rm EXP}}$$  $$u_{{\rm max} \left[1-\left(1+\frac{4t}{\tau _R}\right)\exp \left(-\frac{4t}{\tau _R}\right)\right]H\left(t\right)}$$  $$u_S^{{\rm TNH}}$$  $$\frac{u_{\rm max}}{2} \left[1+\tanh \left(\frac{4t}{\tau _R}\right)\right]$$  STF  Function  $$u_S^{{\rm EXP}}$$  $$u_{{\rm max} \left[1-\left(1+\frac{4t}{\tau _R}\right)\exp \left(-\frac{4t}{\tau _R}\right)\right]H\left(t\right)}$$  $$u_S^{{\rm TNH}}$$  $$\frac{u_{\rm max}}{2} \left[1+\tanh \left(\frac{4t}{\tau _R}\right)\right]$$  View Large The outcome of the parametric analysis is presented herein in terms of comparison with a classical Ground Motion Prediction Equation (GMPE), proposed by Fukushima (2007). The choice of this GMPE is justified by the coherency between the predicted values and aftershock recordings in the Niigata region, although for small magnitude events GMPEs are less controlled by the rupture process, being related more to stress drop and wave propagation effects. Fig. 9 portrays the mentioned comparison in terms of geometric mean of the horizontal PGAs and PGVs. Specifically, Figs 9(a)–(c) pertain the comparison between recordings and the Fukushima GMPEs for PGA and PGV respectively, whereas Figs 9(b)–(d) refer to the SEM3D synthetics. Figure 9. View largeDownload slide Recorded (a) and synthetic (b) horizontal geometric mean of Peak Ground Acceleration (PGAH) obtained by SEM3D simulations of the MJMA4.4 aftershock event, compared to the GMPE proposed by Fukushima (2007). Recorded (c) and synthetic (d) horizontal geometric mean of Peak Ground Velocity (PGVH) obtained by SEM3D simulations of the MJMA4.4 aftershock event, compared to the GMPE proposed by Fukushima (2007). Red circles and blue triangles refer to Brune (uEXP) and a Bouchon’s (uTNH) functions respectively, with free-parameter τR1 = 0.113s; green triangles refer to a Bouchon’s ramp function with free-parameter τR2 = 0.820s instead. Thick black line represents the median PGAH/PGVH stemming from the chosen GMPE; dashed black lines (±σ) represent the standard deviations. Figure 9. View largeDownload slide Recorded (a) and synthetic (b) horizontal geometric mean of Peak Ground Acceleration (PGAH) obtained by SEM3D simulations of the MJMA4.4 aftershock event, compared to the GMPE proposed by Fukushima (2007). Recorded (c) and synthetic (d) horizontal geometric mean of Peak Ground Velocity (PGVH) obtained by SEM3D simulations of the MJMA4.4 aftershock event, compared to the GMPE proposed by Fukushima (2007). Red circles and blue triangles refer to Brune (uEXP) and a Bouchon’s (uTNH) functions respectively, with free-parameter τR1 = 0.113s; green triangles refer to a Bouchon’s ramp function with free-parameter τR2 = 0.820s instead. Thick black line represents the median PGAH/PGVH stemming from the chosen GMPE; dashed black lines (±σ) represent the standard deviations. The Fukushima (2007) GMPE predicts rather well both the recorded mean PGAs and PGVs (±σ), although a slight overestimation may be noticed at larger source-to-site distances. This discrepancy might be due to the magnitude (greater or equal to 5.0) to calibrate the GMPE parameters (Cotton et al.2008). However, Figs 9(b)–(d) suggest that the choice of a rise-time τR1 (for either $$u_S^{{\rm EXP}}$$ or $$u_S^{{\rm TNH}}$$) is the most consistent with the selected GMPE (Fukushima 2007), so it has been chosen hereafter as the suitable STF for the aftershock analyses. It is worth noticing that the choice of a constant Q factor might affect the peak values, diminishing their values compared to the attenuation relationship. 5 SIMULATION OF NEAR-SOURCE WAVEFIELD: SMALL MODEL At this point, the skeleton geology (LARGE model) was upgraded by obtaining an hybrid one (the 3-D structure referred hereafter as SMALL) so to insert the folded structure of the Earth’s crust underneath KKNPP (see the geographical map in Fig. 10a). The latter 3-D geology is defined by seven strata: a representative cross-section SC (see Fig. 10d, corresponding to section S4 in Tsuda et al.2011), passing through KKNPP, was chosen and extruded up to SW (south-westward) and SE (north-eastward). The folding shape is linearly smoothed (the so called smoothing bands in Fig. 10a) so to reconnect with the external subhorizontally layered geology defined in the previous section. This choice was justified by the possible spurious noise generated by the boundaries of the folding zone and polluting the synthetics. Figure 10. View largeDownload slide (a) Geographical map of the KKNPP surroundings, depicting folded cross-section SC, the smoothing bands and the extrusion edged, cross-sections SE and SW. (b and c) Sketched of the smoothed folded structure plugged into the SMALL model: Section SW (b) and Section SC (c) The colours indicate different layers. Figure 10. View largeDownload slide (a) Geographical map of the KKNPP surroundings, depicting folded cross-section SC, the smoothing bands and the extrusion edged, cross-sections SE and SW. (b and c) Sketched of the smoothed folded structure plugged into the SMALL model: Section SW (b) and Section SC (c) The colours indicate different layers. Table 4 summarizes the mechanical properties assumed in this paper for the folding cross-section SC (firstly proposed by Watanabe et al. (2009, 2011); Tsuda et al. (2011)). Table 4. Geological properties of the folding structure underneath KKNPP. VP and VS are the pressure- and shear-wave velocities respectively. The ⋆⋆ indicates the interface chosen to plug the folding structure into the original 1-D Aochi2013 profile, granting a smooth transition from one model to the other. Layer  VS[m s − 1]  VP[m s − 1]  ρ[Kg m − 3]  Nishiyama  700  1900  1700  Shiiya  1200  2200  2100  Upper Teradomari  1700  3300  2300  Lower Teradomari  2000  4200  2400  Nanatani  2000  4600  2500  ⋆⋆Green tuff  2600  5200  2600  Seismic bed rock  2600->2750  5200  2600  Layer  VS[m s − 1]  VP[m s − 1]  ρ[Kg m − 3]  Nishiyama  700  1900  1700  Shiiya  1200  2200  2100  Upper Teradomari  1700  3300  2300  Lower Teradomari  2000  4200  2400  Nanatani  2000  4600  2500  ⋆⋆Green tuff  2600  5200  2600  Seismic bed rock  2600->2750  5200  2600  View Large Although the mechanical properties of the folding model do not differ much from the Aochi2013 profile below 2 km deep, the former geological profile is featured by thinner strata and by lower minimum shear-wave velocity. However, a slight correction in the values of the seismic bedrock shear wave velocity VS was performed herein. The latter was adjusted to grant a smooth transition from the Aochi2013 horizontally layered regional structure to the folding detail. To naturally plug the latter into the former, two welding surfaces were selected, being respectively the upper surface of Aochi2013 layer located at 5 km depth (whose estimated VS value is 2750.0 m s−1, indicated with ⋆ in Table 2) and the lower interface of the Green tuff layer (whose estimated VS value is 2600.0 m s−1, indicated with ⋆⋆ in Table 4) as counterpart in the folding structure. In between the two surfaces, a linear variation of the VS, VP and ρ values is considered. Just replacing the local folding model within the broad Aochi2013 model would produce (as shown by Tsuda et al. (2011) artificial refracted waves, generated around the boundaries/gaps into the hybrid velocity structure. In the SMALL model, moreover, the Japan sea was considered. Due to the asymmetric nature of the folded structural geology beneath the KKNPP site (see Fig. 5a), two small aftershocks scenarios were simulated (see Table 5), whose hypocentres were located along the directrix of the Madonosaka syncline, at the two opposite sides with respect to its planar cross-section passing through the TEPCO facility (see Fig. 10a). Table 5. Summary of the aftershock parameters employed in this analysis. (ϕS; λ; δ) represents the strike, rake and dip angles estimated by F-NET Centroid Moment Tensor solution (Kubo et al.2002). Event  MJMA  M0[Nm]  (ϕS; λ; δ)[°]  τR[s]  AS1 (07/16/07-21:08)  4.4  52.1 · 1014  187; 70; 54  0.113  AS2 (07/16/07-17:42)  4.2  2.09 · 1014  309; 78; 37  0.045  Event  MJMA  M0[Nm]  (ϕS; λ; δ)[°]  τR[s]  AS1 (07/16/07-21:08)  4.4  52.1 · 1014  187; 70; 54  0.113  AS2 (07/16/07-17:42)  4.2  2.09 · 1014  309; 78; 37  0.045  View Large The aim is to characterize the 3-D effect of such a complex geology, by spanning different incidence angles of the wave front impinging the KKNPP. With respect to the LARGE model, the SMALL one, featured by a minimal shear-wave velocity of 700.0 m s−1 was designed to accurately propagate up to 5.0 Hz, with a third order polynomial degree. Except the eventual poor reproduction of the late surface wave arrivals, general features of observed records for AS1 and AS2 could be reproduced, as portrayed in Figs 11 (AS1, 1G1) and 13 (AS2, 5G1). Unfortunately, no recordings for AS1 at 5G1 are in the authors’ disposal. The improvement granted by the inclusion of the folded configuration however is evident (red traces in Figs 11b and d and 13b and d), compared to the poor fit to the records obtained with a subhorizontally layered geology not including the folding (green traces in Figs 11a and c, 12a and c and 13a and c). Figure 11. View largeDownload slide AS1: recorded (REC, blue) and synthetic velocigrams for subhorizontally layered not including the syncline-anticline structures (LAY, green) and folded (FLD, red) geology along the EW (top row) and NS (bottom row) directions, at 1G1. Figure 11. View largeDownload slide AS1: recorded (REC, blue) and synthetic velocigrams for subhorizontally layered not including the syncline-anticline structures (LAY, green) and folded (FLD, red) geology along the EW (top row) and NS (bottom row) directions, at 1G1. Figure 12. View largeDownload slide AS2: recorded (REC, blue) and synthetic velocigrams for subhorizontally layered not including the syncline-anticline structures (LAY, green) and folded (FLD, red) geology along the EW (top row) and NS (bottom row) directions, at 1G1. Figure 12. View largeDownload slide AS2: recorded (REC, blue) and synthetic velocigrams for subhorizontally layered not including the syncline-anticline structures (LAY, green) and folded (FLD, red) geology along the EW (top row) and NS (bottom row) directions, at 1G1. Figure 13. View largeDownload slide AS2: recorded (REC, blue) and synthetic velocigrams for subhorizontally layered not including the syncline-anticline structures (LAY, green) and folded (FLD, red) geology along the EW (top row) and NS (bottom row) directions, at 5G1. Figure 13. View largeDownload slide AS2: recorded (REC, blue) and synthetic velocigrams for subhorizontally layered not including the syncline-anticline structures (LAY, green) and folded (FLD, red) geology along the EW (top row) and NS (bottom row) directions, at 5G1. Figure 14. View largeDownload slide Scores of the Anderon’s Criteria along the EW/NS directions for AS1 (a and b) and AS2 (c and d). Figure 14. View largeDownload slide Scores of the Anderon’s Criteria along the EW/NS directions for AS1 (a and b) and AS2 (c and d). Figure 15. View largeDownload slide AS1: Comparison between recorded (REC, red for EW and orange for NS) response spectra Sa for SG4 and 1G1 sites compared to synthetic ones. (a and c) Sa spectra for layered geology not including the syncline-anticline structure; (b and d) Sa spectra for folded geology. Figure 15. View largeDownload slide AS1: Comparison between recorded (REC, red for EW and orange for NS) response spectra Sa for SG4 and 1G1 sites compared to synthetic ones. (a and c) Sa spectra for layered geology not including the syncline-anticline structure; (b and d) Sa spectra for folded geology. Figure 16. View largeDownload slide AS2: comparison between recorded (REC, red for EW and orange for NS) response spectra Sa for 5G1 and 1G1 sites compared to synthetic ones. (a and c) Sa spectra for layered geology not including the syncline-anticline structure; (b and d) Sa spectra for folded geology. Figure 16. View largeDownload slide AS2: comparison between recorded (REC, red for EW and orange for NS) response spectra Sa for 5G1 and 1G1 sites compared to synthetic ones. (a and c) Sa spectra for layered geology not including the syncline-anticline structure; (b and d) Sa spectra for folded geology. For both AS1 and AS2 and for both Unit 1 and Unit 5, the synthetic waveforms obtained with layered geology (LAY) at surface look out of phase with respect to the records, on one of the two directions EW/NS. This stress the great achievement obtained by including the Madonosaka syncline in the analysis. A more equitable judgment of the waveform fit is provided by the Anderson’s Criteria (Anderson 2004) in Figs 14(a) and (b) for AS1 at 1G1, and in Figs 14(c) and (d) for AS2 at 5G1. For both AS1 and AS2, the complex geological structure seems to influence principally the EW direction. However, the poor score obtained for both cases for the first four Anderson’s Criteria (related to the energy of the waveform) highlights one major drawback of the model: recordings are more energetic than synthetics, due to late arrivals that the model is not fit to reproduce. In terms of normalized Sa response spectra (with respect to their respective PGA), a suitable measure of the goodness of the model is the difference between the natural periods corresponding to the maximum Sa ordinate and referred as to TP. In Figs 15(b)–(d), the focusing of the wavefield due to the syncline causes a shift of TP towards the observed one, at both SG4 (Service Hall, G.L. −250 m) and 1G1 (Unit 1, G.L. 0 m), compared to the poor fit obtained for a subhorizontally layered geology not including the folding (see Figs 15a–c). For AS2, the pseudospectral acceleration response is less sensitive at Unit 5 (Figs 16a and b), whereas the fit is effectively improved at 1G1, by the presence of the folding in the model (Figs 16c and d). Finally, the KKNPP site response estimated by SEM3-D analyses for AS1 and AS2 is condensed in the following four figures: Figs 17 and 18 portray the pseudospectral acceleration spectra Sa estimated at KKNPP, for AS1 and AS2 respectively. Figure 17. View largeDownload slide AS1: Sa response spectra (in cm s−2) at different location around the KKNPP site. LAYERED (red) and FOLDED geological models were compared. Figure 17. View largeDownload slide AS1: Sa response spectra (in cm s−2) at different location around the KKNPP site. LAYERED (red) and FOLDED geological models were compared. Hayakawa et al. (2011) showed that the large amplitude of Unit 1 (located on the southwestward part of the KKNPP) comes from the folding structure where Unit 1 is located on the synclinal axis (Fig. 5a). In agreement with what proposed in the literature, Figs 17 and 18 indicate the folded geology as responsible of higher peaks for the response spectra obtained nearby Unit 1, for natural periods T < 0.5s. The site response is seemingly unaltered by the introduction of a complex folding structure beneath KKNPP (at Unit 5 and Service Hall). The minimum shear velocity introduced in the model corresponds to the engineering bedrock related to depth, justifying that the spatial incoherence of the ground motion recorded at surface and specifically the amplification southwestward are extremely influenced by the syncline–anticline structure, despite the effects due to shallow borehole geology (dispersion, attenuation/amplification) which have been disregarded in the analysis. This strengthens the argument that the site response at the Service Hall (which is located above Chuo-Yutai anticline) depended mainly on the nonlinear site-effects taking place at shallow depths (i.e. G.L.−250 m). The amplification trend occurs independently on the source position, being more accentuated along the EW direction. Watanabe et al. (2009) showed that the Upper Teradomari stratum does not alter the wave propagation of up-going waves, which tend to focalize at the Madonosaka syncline passing through Shiiya stratum. Site amplification becomes therefore significant at Unit 1 (KK1). On the other hand, Unit 5 (KK5, which is located above the Ushirodani anticline) is evidently more sensitive to wave-motion travelling from southwest to northeast, throughout the folding zone: the pseudospectral peaks differ from layered to folded geology only in Fig. 18, referring to AS2 (nucleated close by the third asperity on the fault rupture plane). Figure 18. View largeDownload slide Sa response spectra (in cm s−2) at different location around the KKNPP site. LAYERED (red) and FOLDED geological models were compared. Figure 18. View largeDownload slide Sa response spectra (in cm s−2) at different location around the KKNPP site. LAYERED (red) and FOLDED geological models were compared. 6 ON THE ACCURACY OF NUMERICAL SOLUTION The accuracy of the numerical method employed to solve the 3-D elastodynamic problem is notoriously a big deal to cope with. Due to inherent spatiotemporal nature of the wave-propagation physical phenomenon, the origins of the observed numerical dispersion are namely (1) the spatial and the (2) time discretization indeed. Usually, a Finite-Difference of Finite-Element spatial semi-discretization is assumed, along with a numerical time-integrator belonging to the Newmark’s family. In the SEM, high-order polynomials are employed, sampled at the GLL point on each element. Despite the intrinsic h − p property held by the SEM (leading to finer spatial discretizations and increased accuracy), it is intuitive that the two sources of numerical dispersions are not disconnected: Seriani & Oliveira (2008) stated that low-order time discretizations (such as the 2nd-order accurate Newmark’s schemes) deteriorate the high accuracy in space. According to Seriani & Oliveira (2008), P-wave dispersion and polarization errors are less sensitive to the value of Poisson’s ratio than the S-wave dispersion error. They confirmed that four grid points per wavelength are sufficient to have the dispersion error below 1 per cent on SE approximations of degree eight with GLL collocation points. Following the indications provided by the modal analysis performed by the mentioned authors, a parametric analysis on the computational model presented in this study for folded geology was performed. The minimum element size of the computational grid is ≈250 m at surface, with VS, min = 700.0 m s−1 and Poisson’s ratio ν ≈ 0.4. Three analysis were run therefore, with 5, 7 and 10 GLL points per element edge respectively. According to Seriani & Oliveira (2008) and for the Poisson’s ratio considered, the frequencies delimiting the accuracy of the analyses are 2.8, 5.0 and 7.0 Hz respectively. Therefore, the whole set of results of the three-analyses were low-pass filtered at 7 Hz, to check whether significant differences emerge. Fig. 19 shows the comparisons (in terms of accelerograms, Fourier’s spectrum along NS and EW directions respectively and geometric mean of the horizontal pseudospectral acceleration SaGMH) for 5, 7 and 10 GLL points respectively, at 1G1 (Unit 1, G.L. 0 m). The comparisons are rather explicative, in the following sense: the ground motion estimation does not change significantly for 7 or 10 GLL points, whereas it is rather underestimated (see the SaGMH ordinates in Fig. 19c) when the computational model is featured by 5 GLL points; the response at short periods is amplified by increasing the number of GLL points, granting a certain confidence to the analyses portrayed in the previous section (comparison with layered geology). Increasing the number of GLL points would lead to even more amplified response at both 1G1 and 5G1. the computational costs (measured in CPU-time, and plotted in the captioned subaxes in Fig. 19c) are growing exponentially with increasing number of GLL points, although no remarkable difference can be found between 7 and 10 GLL. This justify, for this particular case, the choice of 7 GLL as an optimum solution for the range of frequency of interest. Figure 19. View largeDownload slide Results of the parametric analysis on the SMALL model with folded geology at 1G1. (a) Accelerograms(NS direction), (b) Fourier’s spectra (NS direction) and (c) pseudospectral acceleration responses (geometric mean of horizontal components) for 5 (blue), 7 (red) and 10 (green) GLL respectively. Synthetics were filtered at 7 Hz. The computational costs of the SEM3D analyses are plotted in the captioned axes in (c). Figure 19. View largeDownload slide Results of the parametric analysis on the SMALL model with folded geology at 1G1. (a) Accelerograms(NS direction), (b) Fourier’s spectra (NS direction) and (c) pseudospectral acceleration responses (geometric mean of horizontal components) for 5 (blue), 7 (red) and 10 (green) GLL respectively. Synthetics were filtered at 7 Hz. The computational costs of the SEM3D analyses are plotted in the captioned axes in (c). 7 CONCLUSIONS In this paper, a source-to-site numerical model of the NCOEQ-2007 earthquake is constructed, to reproduce the seismic site response observed at KKNPP. Specifically, the aim of this study is to criticize the common simplified assumption of a subhorizontally layered Earth’s crust, inferred from seismic and geodetic and geological surveys, stressing the importance of the shallow geological conformation in altering the incident wave motion. The case of KKNPP rises as an exemplary benchmark, due to extensively observed near-field effects at the site, such as the ground motion incoherence within a relatively small distance. As a matter of fact, the incident wave motion resulted amplified in the south-west area, nearby the first group of nuclear reactors and structures (Unit 1–4), compared to the second group of facilities (Unit 5–7), located northeastward. Two aftershocks belonging to the NCOEQ-2007 sequence were simulated to prove the assumption focusing on the effect of the geological configuration. The dense observation network deployed in the Niigata region and at the KKNPP steered the calibration of two numerical seismic scenarios: one to portrays the regional incident wavefield (in a frequency band 0.1–0.5 Hz) and a second one restricted to the ground shaking pertaining the KKNPP surroundings, but extending the modelling accuracy up to 5.0 Hz (and possibly 7.0 Hz). The paper explains the methodological approach to construct an holistic source-to-site computational model: LARGE model is first calibrated upon the results of semi-analytical solutions, considering a 1-D regional geology so to constrain the low-frequency band of the incident wavefield; SMALL model is later constructed by trimming a portion of LARGE model and plugging into the 1-D regional geology the intricate 3-D set of syncline-anticlines, placed right underneath KKNPP. LARGE model is used to assess the effect of the source parameters and of the topography, comparing the results with classical GMPEs. SMALL model includes the Japan Sea, the bathymetry and the coastline, other than the folding geology. The outcome of the numerical exercise provides a fairly good agreement between the recordings and the synthetics at KKNPP, with high accuracy at 5.0 Hz. The effect of the folding clearly improves this fit, compared to the too much simplified subhorizontally layered profiles. Moreover, the effort of including a complicated geology in the seismological model is worthwhile, as already proved by Tsuda et al. (2011), but herein presented in a more extended way, including rigorous comparison between time-histories, goodness of fit estimation at several locations and detailed seismic site response at KKNPP. One major contribution of this paper is that the two-step construction of a hybrid velocity structure bridges the gap between previous studies targeting the regional wavefield (but limited in frequency, e.g. Aochi et al.2013) and more site-specific analyses, describing the wavefield at KKNPP (Watanabe et al.2009; Tsuda et al.2011). The numerical analyses highlight the impact of the syncline laying below Unit 1 in causing an amplification (both in terms of peak ground motions and response spectra) with respect to the layered geology. The syncline seemingly focalizes the upgoing wavefield, drifting the radiated energy towards Unit 1. This occurrence represents an exception to the current design standards, which traditionally leans towards the assumption of a 1-D layered geology in case of lack of further information concerning the subsurface geology. The case study stresses the importance of the propagation path in drifting and altering (eventually amplifying) the wavefield radiated even from relatively weak aftershocks. Those findings confirm previous studies and extend their validity at higher frequency (5.0 Hz). Moreover, in this paper a major effort of quantifying the numerical accuracy is performed and presented aside the main analysis, completing the results and depicting the uncertainty related. However, the major shortcoming of the analyses presented herein resides in the lack of more refined models of the scattered wavefield in a high frequency range (0–10 Hz). This is a crucial feature to be taken into consideration in the transition towards a broad-band deterministic modelling of the earthquake phenomenon. For instance, the inclusion of heterogeneous soil deposits and nonlinear hysteretic damping would surely enhance the realism of the final outcome, although it may mislead the understanding of the physical mechanisms taking place, due to the increased number of parameters required. Based on these considerations, it can be argued that the current seismic design guidelines for critical structures should require further improvements of site-characterization not only in terms of active faults and local geotechnical conditions, but specifically on the conformation of the Earth’s crust (such as the SCEC research plan for the period 2012-2016, focused on specific areas in the Southern California and enhanced the Community Velocity Model (CVM), describing seismic P- and S-wave velocities and densities throughout the southern California region). It is in the authors’ belief that the modern procedures for earthquake ground motion prediction for seismic hazard analysis, especially for the low probability events inherent to the seismic design of critical structures, should account for a reliable incident wavefield issued from forward deterministic analyses. The ever increasing computational power leads to put major efforts to extend the frequency band and to enable PBS to be used with confidence in engineering applications (e.g. seismic hazard analysis and structural or geotechnical dynamic analysis). Acknowledgements This work, within the SINAPS@ project, benefited from French state funding managed by the National Research Agency under program RNSR Future Investments bearing reference No. ANR-11-RSNR-0022-04. The research reported in this paper has been supported in part by the SEISM Paris Saclay Research Institute. The authors are very grateful to the Tokyo Electric Power Company (TEPCO) for providing such high-quality earthquake recordings of the Niigata-Ken Chūetsu-Oki, reported in TEPCO (2007). Time histories and velocity profiles used in this study were collected from the KiK-net website: http://www.kik.bosai.go.jp/kik/ (last accessed November 2011). The authors are very grateful to the National Research Institute for Earth Science and Disaster Prevention (NIED) for providing such high-quality earthquake recordings. The authors are very thankful to the Data Management Center of the National Research Institute for Earth Science and Disaster Resilience (NIED DMC), which provides continuous seismic waveform and other seismological data obtained by the NIED High Sensitivity Seismograph Network (NIED Hi-net) and the other seismograph networks in Japan. NIED, Japan Meteorological Agency (JMA), and universities have started exchanging high-sensitivity seismograph data obtained by each seismic network with a specific policy for ‘Exchange and Open Access of Waveform Data obtained by High Sensitivity Seismograph Observation’ conducted by the Headquarters for Earthquake Research Promotion of the Japanese government. The website (http://www.hinet.bosai.go.jp/?LANG = en) provides not only the NIED Hi-net waveform data but also the other high-sensitivity seismograph data. The Digital Elevation Model belongs to the Shuttle Radar Topography Mission (SRTM) database, result from a collaborative effort by the National Aeronautics and Space Administration (NASA) and the National Geospatial-Intelligence Agency (NGA), as well as the participation of the German and Italian space agencies, to generate a near-global digital elevation model (DEM) of the Earth using radar interferometry. Source: https://dds.cr.usgs.gov/srtm/version2_1. REFERENCES Aki K., 1967. Scaling law of seismic spectrum, J. geophys. Res. , 72( 4), 1217– 1231. https://doi.org/10.1029/JZ072i004p01217 Google Scholar CrossRef Search ADS   Anderson J.G., 2004. Quantitative measure of the goodness-of-fit of synthetic seismograms, in13th World Conference on Earthquake Engineering, Vancouver, BC, Canada, Paper No. 243, pp. 775– 784. Aochi H., Yoshimi M., 2016. Seismological asperities from the point of view of dynamic rupture modeling: the 2007 Mw6.6 Chuetsu-Oki, Japan, earthquake, J. Seismol. , 20( 4), 1089– 1105. https://doi.org/10.1007/s10950-016-9569-5 Google Scholar CrossRef Search ADS PubMed  Aochi H., Ducellier A., Dupros F., Delatre M., Ulrich T., de Martin F., Yoshimi M., 2013. Finite difference simulations of seismic wave propagation for the 2007 Mw6.6 Niigata-ken Chuetsu-Oki earthquake: validity of models and reliable input ground motion in the near field, Pure appl. Geophys. , 170( 1–2), 43– 64. https://doi.org/10.1007/s00024-011-0429-5 Google Scholar CrossRef Search ADS   Baker J.W., 2007. Quantitative classification of near-fault ground motions using wavelet analysis, Bull. seism. Soc. Am. , 97( 5), 1486– 1501. https://doi.org/10.1785/0120060255 Google Scholar CrossRef Search ADS   Bao H., Bielak J., Ghattas O., Kallivokas L., O’Hallaron D., Shewchuk J., Xu J., 1998. Large-scale simulation of elastic wave propagation in heterogeneous media on parallel computers, Comput. Methods Appl. Mech. Eng. , 152( 1), 85– 102. https://doi.org/10.1016/S0045-7825(97)00183-7 Google Scholar CrossRef Search ADS   Bard P.-Y., Bouchon M., 1980. The seismic response of sediment-filled valleys. Part 1. The case of incident SH waves, Bull. seism. Soc. Am. , 70( 4), 1263– 1286. Bizzarri A., 2014. On the point-source approximation of earthquake dynamics, Ann. Geophys. , 57( 3). doi: http://dx.doi.org/10.4401/ag-6479. Bouchon M., 1981. A simple method to calculate Green’s functions for elastic layered media, Bull. seism. Soc. Am. , 71( 4), 959– 971. Brune J.N., 1970. Tectonic stress and the spectra of seismic shear waves from earthquakes, J. geophys. Res. , 75( 26), 4997– 5009. https://doi.org/10.1029/JB075i026p04997 Google Scholar CrossRef Search ADS   Chaljub E. et al.  , 2015. 3-D numerical simulations of earthquake ground motion in sedimentary basins: testing accuracy through stringent models, Geophys. J. Int. , 201( 1), 90– 111. https://doi.org/10.1093/gji/ggu472 Google Scholar CrossRef Search ADS   Cirella A., Piatanesi A., Tinti E., Cocco M., 2008. Rupture process of the 2007 Niigataken Chuetsu-Oki earthquake by non-linear joint inversion of strong motion and GPS data, Geophys. Res. Lett. , 35( 16), 1– 5. https://doi.org/10.1029/2008GL034756 Google Scholar CrossRef Search ADS   Cotton F., Pousse G., Bonilla F., Scherbaum F., 2008. On the discrepancy of recent European ground-motion observations and predictions from empirical models: Analysis of KiK-net accelerometric data and point-sources stochastic simulations, Bull. seism. Soc. Am. , 98( 5), 2244– 2261. https://doi.org/10.1785/0120060084 Google Scholar CrossRef Search ADS   Cox K.E., Ashford S.A., 2002. Characterization of large velocity pulses for laboratory testing, Tech. Rep. April , Pacific Earthquake Engineering Research Center. Day S.M., Graves R., Bielak J., Dreger D., Larsen S., Olsen K.B., Pitarka A., Ramirez-Guzman L., 2008. Model for basin effects on long-period response spectra in southern California, Earthq. Spectra , 24( 1), 257– 277. https://doi.org/10.1193/1.2857545 Google Scholar CrossRef Search ADS   De Martin F., 2011. Verification of a spectral-element method code for the Southern California earthquake center LOH.3 viscoelastic case, Bull. seism. Soc. Am. , 101( 6), 2855– 2865. https://doi.org/10.1785/0120100305 Google Scholar CrossRef Search ADS   De Martin F., Aochi H., Foerster E., 2007. Testing the Double-Couple Point Source Model in GEFDYN, Tech. rep. , BRGM. Dreger D., Hurtado G., Chopra A.K., Larsen S., 2007. Near-Fault Seismic Ground Motions, Tech. Rep. 59A0435 , California Department of Transportation. Ducellier A., Aochi H., 2010. Numerical simulation of the Mw6.6 Niigata, Japan, earthqhake: reliable input ground motion for engineering purpose, in 14th European Conference on Earthquake Engineering , Ohrid, Macedonia. Dupros F., De Martin F., Foerster E., Komatitsch D., Roman J., 2010. High-performance finite-element simulations of seismic wave propagation in three-dimensional nonlinear inelastic geological media, Parallel Comput. , 36( 5-6), 308– 325. https://doi.org/10.1016/j.parco.2009.12.011 Google Scholar CrossRef Search ADS   Duputel Z., Rivera L., Kanamori H., Hayes G., 2012. W phase source inversion for moderate to large earthquakes (1990-2010), Geophys. J. Int. , 189( 2), 1125. https://doi.org/10.1111/j.1365-246X.2012.05419.x Google Scholar CrossRef Search ADS   Duputel Z., Tsai V.C., Rivera L., Kanamori H., 2013. Using centroid time-delays to characterize source durations and identify earthquakes with unique characteristics, Earth planet. Sci. Lett. , 374, 92– 100. https://doi.org/10.1016/j.epsl.2013.05.024 Google Scholar CrossRef Search ADS   Faccioli E., Vanini M., 2003. Complex seismic site effects in sediment-filled valleys and implications on design spectra, Prog. Struct. Eng. Mater. , 5( 4), 223– 238. https://doi.org/10.1002/pse.156 Google Scholar CrossRef Search ADS   Faccioli E., Maggio F., Paolucci R., Quarteroni A., 1997. 2D and 3D elastic wave propagation by a pseudo-spectral domain decomposition method, J. Seismol. , 1( 3), 237– 251. https://doi.org/10.1023/A:1009758820546 Google Scholar CrossRef Search ADS   Fujiwara H. et al.  , 2009. A study on subsurface structure model for deep sedimentary layers of Japan for strong-motion evaluation, Technical Note of the National Research Institute for Earth Science and Disaster Prevention , No. 337, 260. Fukushima S., Hayashi T., Yashiro H., 2007. Seismic hazard analysis based on the joint probability density function of PGA and PGV, in Transactions, SMiRT 19 , Vol. M03. Gallovič F., Käser M., Burjánek J., Papaioannou C., 2010. Three-dimensional modeling of near-fault ground motions with nonplanar rupture models and topography: case of the 2004 Parkfield earthquake, J. geophys. Res. , 115( B3), 1– 17. Google Scholar CrossRef Search ADS   Gélis C., Bonilla L.F., 2012. 2-D P-SV numerical study of soil-source interaction in a non-linear basin, Geophys. J. Int. , 191, 1374– 1390. Göddeke D., Komatitsch D., Möller M., 2014. Finite and spectral element methods on unstructured grids for flow and wave propagation methods, chap.9, in Numerical Computations with GPUs , pp. 183– 206, ed. Kindratenko V., Springer. Graves R.W., 1996. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences, Bull. seism. Soc. Am. , 86( 4), 1091– 1106. Graves R.W., Wald D.J., 2004. Observed and simulated ground motions in the San Bernardino Basin Region for the Hector Mine, California, earthquake, Bull. seism. Soc. Am. , 94( 1), 131– 146. https://doi.org/10.1785/0120030025 Google Scholar CrossRef Search ADS   Gürpinar A., Serva L., Livio F., Rizzo P.C., 2017. Earthquake-induced crustal deformation and consequences for fault displacement hazard analysis of nuclear power plants, Nucl. Eng. Des. , 311, 69– 85. https://doi.org/10.1016/j.nucengdes.2016.11.007 Google Scholar CrossRef Search ADS   Hartzell S.H., 1978. Earthquake aftershocks as Green’s functions, Geophys. Res. Lett. , 5( 1), 1– 4. https://doi.org/10.1029/GL005i001p00001 Google Scholar CrossRef Search ADS   Hayakawa T., Tsuda K., Uetake T., Hikima K., Tokumitsu R., Nagumo H., 2011. Modeling 3D velocity structure in the fault region of the 2007 Niigataken Chuetu-oki Earthquake-Incorporating the 3D fold geological structure beneath the Kashiwazaki-Kariwa nuclear power plant, in Japan Geoscience Union Meeting , no. SSS023, p. 14. Hijikata K., Nishimura I.and Mizutani H., Tokumitsu R., Mashimo M., Tanaka S., 2010. Ground motion characteristics of 2007 Niigata-ken Chuetsu-Oki earthquake, J. Struct. Construct. Eng. , 75( 653), 1279– 1288. https://doi.org/10.3130/aijs.75.1279 Google Scholar CrossRef Search ADS   Hisada Y., 1994. An efficient method for computing Green’s Functions for layered half-space with sources and receivers at close depths, Bull. seism. Soc. Am. , 84( 5), 1456– 1472. Hisada Y., 1995. An efficient method for computing Green’s functions for a layered half-space with sources and receivers at close depths (part 2), Bull. seism. Soc. Am. , 85( 4), 1080– 1093. Hisada Y., 2008. Broadband strong motion simulation in layered half-space using stochastic Green’s function technique, J. Seismol. , 12( 2), 265– 279. https://doi.org/10.1007/s10950-008-9090-6 Google Scholar CrossRef Search ADS   Irikura K., 1983. Semi-Empirical Estimation of Strong Ground Motions During, Bull. Disaster Prevention Res. Inst. , 33( 2), 63– 104. JNES, 2008. Construction of the subsurface structure around the source area of 2007 Niigata Ken Chuetsu-Oki earthquake, Tech. rep. , Japanes Nuclear Energy Safety Organization. Kato A., Sakai S., Kurashimo E., Igarashi T., Iidaka T., Hirata N., Iwasaki T., Kanazawa T., 2008. Imaging heterogeneous velocity structures and complex aftershock distributions in the source region of the 2007 Niigataken Chuetsu-Oki Earthquake by a dense seismic observation, Earth Planets Space , 60, 1111– 1116. https://doi.org/10.1186/BF03353145 Google Scholar CrossRef Search ADS   Kato A. et al.  , 2009. Reactivation of ancient rift systems triggers devastating intraplate earthquakes, Geophys. Res. Lett. , 36( L05301), 1– 5. Kawabe H., Kamae K., 2010. Source modelling and 3D ground motion simulation of the 2007 Niigataken Chuetsu-oki earthquake (Mj6.8), in 13th Japan Earthquake Engineering Symposium , Tsukuba Japan. Kayen R. et al.  , 2007. Investigation of the M6.6 Niigata-Chuetsu Oki, Japan, Earthquake of July 16, 2007, Tech. Rep. 1365 , U.S. Department of the Interior and U.S. Geological Survey (USGS). Kobayashi I., Tateishi M., Yoshimura T., Ueda T., Kato H., 1995. Geology of the Kashiwazakii district. Quadrangle Series: Scale 1:50000, Vol. 7, Issue 37, P. 101, Geological Survey of Japan. Komatitsch D., Vilotte J.-P., 1998. The spectral element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures, Bull. seism. Soc. Am. , 88( 2), 368– 392. Komatitsch D., Vilotte J.P., Vai R., Castillo-Covarrubias J.M., Sánchez-Sesma F.J., 1999. The spectral element method for elastic wave equations: application to 2D and 3D seismic problems, in SEG Expanded Abstracts 17 , vol. 45, pp. 1139– 1164. Kubo A., Fukuyama E., Kawai H., Nonomura K., 2002. NIED seismic moment tensor catalogue for regional earthquakes around Japan: quality test and application, Tectonophysics , 356( 1), 23– 48. https://doi.org/10.1016/S0040-1951(02)00375-X Google Scholar CrossRef Search ADS   Mazzieri I., Smerzini C., Antonietti P.F., Rapetti F., Stupazzini M., Paolucci R., Quarteroni A., 2011. Non-conforming spectral approximations for the elastic wave equation in heterogeneous media, in COMPDYN 2011 ECCOMAS Thematic Conference on Computational Methods in Structural Dynamics and Earthquake Engineering , pp. 1– 17, eds Papadrakakis M., Fragiadakis M., Plevris V., Corfu, Greece. Mazzieri I., Stupazzini M., Guidotti R., Smerzini C., 2013. SPEED: SPectral Elements in Elastodynamics with Discontinuous Galerkin: a non-conforming approach for 3D multi-scale problems, Int. J. Numer. Methods Eng. , 95( 12), 991– 1010. https://doi.org/10.1002/nme.4532 Google Scholar CrossRef Search ADS   Miyake H., Koketsu K., Hikima K., Shinohara M., Kanazawa T., 2010. Source Fault of the 2007 Chuetsu-Oki, Japan, Earthquake, Bull. seism. Soc. Am. , 100( 1), 384– 391. https://doi.org/10.1785/0120090126 Google Scholar CrossRef Search ADS   Moczo P., Kristek J., Galis M., Pazak P., Balazovjech M., 2007. The finite-difference and finite-element modeling of seismic wave propagation and earthquake motion, Acta Phys. Slovaca , 57( 2), 177– 406. https://doi.org/10.2478/v10155-010-0084-x Google Scholar CrossRef Search ADS   Obara K., Kasahara K., Hori S., Okada Y., 2005. A densely distributed high-sensitivity seismograph network in Japan: Hi-net by National Research Institute for Earth Science and Disaster Prevention, Rev. Sci. Instrum. , 76( 2), 1– 12. https://doi.org/10.1063/1.1854197 Google Scholar CrossRef Search ADS   Paolucci R., Mazzieri I., Smerzini C., Stupazzini M., 2014. Physics -Based Earthquake Ground Shaking Scenarios in Large Urban Areas, in Perspectives on European Earthquake Engineering and Seismology , Vol. 34 of Geotechnical, Geological and Earthquake Engineering, pp. 331– 359, ed. Ansal A., Springer. Google Scholar CrossRef Search ADS   Paolucci R., Mazzieri I., Smerzini C., 2015. Anatomy of strong ground motion: near-source records and 3D physics-based numerical simulations of the Mw 6.0 May 29 2012 Po Plain earthquake, Italy, Geophys. J. Int. , 203, 2001– 2020. https://doi.org/10.1093/gji/ggv405 Google Scholar CrossRef Search ADS   Pavlenko O.V., Irikura K., 2012. Nonlinear soil behavior at the Kashiwazaki-Kariwa Nuclear Power Plant during the Niigata Chuetsu-Oki earthquake (July, 16, 2007), Pure appl. Geophys. , 169( 10), 1777– 1800. https://doi.org/10.1007/s00024-011-0447-3 Google Scholar CrossRef Search ADS   Pitarka A., Irikura K., Iwata T., Sekiguchi H., 1998. Three-dimensional simulation of the near-fault ground motion for the 1995 Hyogo-Ken Nanbu (Kobe), Japan, earthquake, Bull. seism. Soc. Am. , 88( 2), 428– 440. Quinay P. E.B., Ichimura T., Hori M., Nishida A., Yoshimura S., 2013. Seismic Structural Response Estimates of a Fault-Structure System Model with Fine Resolution Using Multiscale Analysis with Parallel Simulation of Seismic-Wave Propagation, Bull. seism. Soc. Am. , 103( 3), 2094– 2110. https://doi.org/10.1785/0120120216 Google Scholar CrossRef Search ADS   Sekiguchi H. et al.  , 2009. 3D subsurface structure model of the Niigata sedimentary basin, Ann. Rep. Active Fault Paleoearthq. Res. , 9, 175– 259. Seriani G., Oliveira S.P., 2008. Dispersion analysis of spectral element methods for elastic wave propagation, Wave Motion , 45( 6), 729– 744. https://doi.org/10.1016/j.wavemoti.2007.11.007 Google Scholar CrossRef Search ADS   Shiba Y., 2008. Source process and broadband strong motions during the Niigata-ken Chuetsu-Oki earthquake in 2007, Denryoku Chuo Kenkyusho Hokoku , ( N08007), 1– 4, Japanese with English Abstract. Shinohara M. et al.  , 2008. Precise aftershock distribution of the 2007 Chuetsu-Oki Earthquake obtained by using an ocean bottom seismometer network, Earth Planets Space , 60( 11), 1121– 1126. https://doi.org/10.1186/BF03353147 Google Scholar CrossRef Search ADS   Smerzini C., Villani M., 2012. Broadband Numerical Simulations in Complex Near-Field Geological Configurations: The Case of the 2009 Mw 6.3 L’Aquila Earthquake, Bull. seism. Soc. Am. , 102( 6), 2436– 2451. https://doi.org/10.1785/0120120002 Google Scholar CrossRef Search ADS   Spudich P., Archuleta R., 1987. Techniques for earthquake ground-motion calculation with applications to source parameterization to finite faults, in Seismic Strong Motion Synthesis , chap. pp. 205– 265, ed. Bolt B.A., Academic Press. Google Scholar CrossRef Search ADS   Taborda R., Bielak J., 2011. Large-scale earthquake simulation: computational seismology and complex engineering systems, Comput. Sci. Eng. , 13( 4), 14– 27. https://doi.org/10.1109/MCSE.2011.19 Google Scholar CrossRef Search ADS   Taborda R., Bielak J., 2013. Ground-motion simulation and validation of the 2008 Chino Hills, California, earthquake, Bull. seism. Soc. Am. , 103( 1), 131– 156. https://doi.org/10.1785/0120110325 Google Scholar CrossRef Search ADS   TEPCO, 2007. The data analysis recorded at the Kashiwazaki Kariwa Nuclear Power Plant during the 2007 Niigata-ken Chuetsu-oki earthquake, Tech. rep. , The Tokyo Electric Power Company, Inc, in Japanese. Tokumitsu R., Kikuchi M., Nishimura I., Shiba Y., Tanaka S., 2009. Analysis of the strong motion records obtained from the 2007 Niigataken Chuetsu-Oki earthquake and determination of the design basis ground motions at the Kashiwazaki Kariwa Nuclear Power Plant. Part 1. Outline of the strong motion records and estimation of factors in large amplification, Tech. rep. Tsuda K., Hayakawa T., Uetake T., Hikima K., Tokimitsu R., Nagumo H., Shiba Y., 2011. Modeling 3D Velocity Structure in the Fault Region of the 2007 Niigataken Chuetu-Oki Earthquake with Folding Structure, in 4th IASPEI/IAEE International Symposium on Effects of Surface Geology on Seismic Motion , Santa Barbara, CA, pp. 1– 11. Uetake T., Nishimura I., Mizutani H., 2008. Characteristics of strong motion records in Kashiwazaki-Kariwa Nuclear Power Station during the Niigataken Chuetsu-Oki earthquake in 2007, in 14th WCEE , Beijing, China. Watanabe T., Moroi T., Nagano M., Tokumitsu R., Kikuchi M., Nishimura I., 2009. Analysis of the strong motion records obtained from the 2007 Niigataken Chuetsu-Oki earthquake and determination of the design basis ground motions at the Kashiwazaki Kariwa Nuclear Power Plant. Part 2. Difference of site amplification based on the 2D FEM analysis of the folded structure, Tech. rep. Watanabe T., Moroi T., Tokumitsu R., Nishimura I., Hijikata K., 2011. Examination of relation between locations of asperities and site amplification characteristics of ground motions by analysis considering the folded structure-Estimation based on the strong motion records obtained from the 2007 Niigataken Chuetsu-oki earthquake in the Kashiwazaki-Kariwa Nuclear Power Station, J. Struct. Construct. Eng. , 76( 659), 71– 78. https://doi.org/10.3130/aijs.76.71 Google Scholar CrossRef Search ADS   Yee E., Stewart J., Tokimatsu K., 2011. Nonlinear Site Response and Seismic Compression at Vertical Array Strongly Shaken by the Niigata-ken Chuetsu-Oki Earthquake, Tech. rep. , Pacific Earthquake Engineering Research Center-College of Engineering, University of California, Berkeley. APPENDIX: To assess the goodness of fit, the Anderson’s criteria are employed (Anderson 2004, listed in Table A1). Table A1. Goodness of fit parameters. Number  Symbol  Similarity of:  C1  SDa  Arias duration  C2  SDe  Energy duration  C3  SIa  Arias Intensity  C4  SIv  Energy Integral  C5  Spga  Peak Acceleration  C6  Spgv  Peak Velocity  C7  Spgd  Peak Displacement  C8  Ssa  Response Spectra  C9  Sfs  Fourier Spectra  C10  C*  Cross Correlation  Number  Symbol  Similarity of:  C1  SDa  Arias duration  C2  SDe  Energy duration  C3  SIa  Arias Intensity  C4  SIv  Energy Integral  C5  Spga  Peak Acceleration  C6  Spgv  Peak Velocity  C7  Spgd  Peak Displacement  C8  Ssa  Response Spectra  C9  Sfs  Fourier Spectra  C10  C*  Cross Correlation  View Large © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# On the effect of the 3-D regional geology on the seismic design of critical structures: the case of the Kashiwazaki-Kariwa Nuclear Power Plant

, Volume 213 (2) – May 1, 2018
20 pages

/lp/ou_press/on-the-effect-of-the-3-d-regional-geology-on-the-seismic-design-of-KDL6Y5nIf0
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggy027
Publisher site
See Article on Publisher Site

### Abstract

Summary In this study, numerical investigation is performed on a realistic source-to-site earthquake scenario, with the aim to assess the role of complex 3-D geological structures on the predicted wavefield. With this respect, the paper pointedly targets the seismic response of nuclear power plants in near-field conditions and the verification of some simplified assumptions commonly adopted for earthquake ground motion prediction and site effects analysis. To this purpose, the Kashiwazaki-Kariwa Nuclear Power Plant (Japan) is assumed as reference case-study. In 2007, the nuclear site and its surroundings were struck by the Niigata-Ken Chūetsu-Oki seismic sequence, which caused some of the peak ground motion design limits to be largely overpassed. The dense observation network deployed at the site recorded a highly incoherent and impulsive earthquake ground motion. Many studies argued that the intricate syncline-anticline geology lying underneath the nuclear facility was highly responsible of the observed seismic response. Therefore, a physics-based numerical model of the epicentral area is built-up (≈60 km wide) and tested for small aftershocks, so to discount the effect of extended source on the synthetic site-response. The numerical model (based on the Spectral Element Method) reproduces the source-to-site wave propagation by embracing the effects of the surface topography along with the presence of the Japan Sea (i.e. the bathymetry, the coastline and the fluid–solid interaction). Broad-band (0–5 Hz) synthetic waveforms are obtained for two different aftershocks, located at the two opposite sides of the nuclear facility, aiming to assess the influence of the incidence angle the radiated wave field impinges the foldings beneath it. The effect of the folding presence is assessed by comparing it to a subhorizontally layered geology, in terms of numerical outcome, and by highlighting the differences with respect to the observations. The presence of an intricate geology effectively unveils the reason behind the observed ground motion spatial variability within a relatively small area, stressing its crucial role to properly reproduce the modification the wavefield undergoes during its propagation path towards the surface. The accuracy of the numerical exercise is discussed along with its results, to show the high-fidelity of these deterministic earthquake ground motion predictions. Computational seismology, Earthquake dynamics, Earthquake ground motions, Wave propagation 1 INTRODUCTION In the past two decades, the seismic hazard analysis and vulnerability assessment of large urban areas and critical structures and infrastructures have progressively taken advantage of the ever-increasing computational power available (Paolucci et al.2014). This outstanding technological and numerical progress seemingly broke through the evergreen and most stringent bottleneck in computational seismology: the impossibility to solve the complete source-to-site seismic wave propagation problem in a single-step analysis. This achievement paved the way to fully couple the large scale seismological models for the region of interest (i.e. the fault mechanism and the geological properties of the Earth’s crust), with local engineering models for geotechnical, site-effect and structural analyses, in the next few years. To this end, a plethora of extremely efficient numerical methods, for instance, the Finite Difference Method (FDM; Graves 1996; Pitarka et al.1998; Moczo et al.2007), the Finite Element Method (FEM; Taborda & Bielak 2011, 2013) and the Spectral Element Method (SEM; Faccioli et al.1997; Komatitsch & Vilotte 1998; Mazzieri et al.2011) has been progressively applied to the seismological problem. In spite of the inherent complexity and the huge dimensions of those computational models, their power is essentially embodied by the higher broad-band accuracy they provide (i.e. up to 4-5 Hz, De Martin 2011; Paolucci et al.2015), gradually bridging the gap between low-frequency source models obtained via waveform inversion techniques and the structural modal frequencies (i.e. up to 20 Hz). The mentioned enhancement in computational seismology was certainly granted by improved vectorization strategies to implement and tailor the mentioned numerical methods to massively parallel architectures (Dupros et al.2010; Mazzieri et al.2013; Göddeke et al.2014). Once shattered the computational barriers, a new holistic philosophy took place, driven by the deterministic modelling of the physics underlying each aspect of the earthquake phenomenon, for more accurate sensitivity analyses and uncertainty quantification of models and related parameters. All the ingredients (i.e. source, path and site effects) are naturally convolved in one-step all-embracing analysis, which is capable to predict a realistic seismic wavefield and to explain the observed time-histories in sedimentary deposits and/or at a global scale (De Martin 2011). A 20 yr history of 3-D physics-based simulations of past earthquakes has been written so far, with an overall increasing frequency upper limit (Graves 1996; Pitarka et al.1998; Bao et al.1998; Komatitsch et al.1999; Graves & Wald 2004; Day et al.2008; Gallovič et al.2010; Smerzini & Villani 2012; Aochi et al.2013; Paolucci et al.2015,among others). For a detailed excursus, please refer to Paolucci et al. (2015), who listed a selection of recent studies to produce earthquake ground shaking scenarios in large metropolitan areas. To clarify, by regional scale, it is hereafter assumed an average distance range between about 1 km (local scale) to 100 km (continental scale). Moreover, several numerical benchmarks have been recently performed during the past years, with the aim at testing the accuracy and the efficiency of different simulation methods when handling complex domain geometries and material interfaces. For instance, Chaljub et al. (2015) tested different 3-D numerical methods for earthquake ground shaking scenarios. They defined four canonical models inspired by the Mygdonian basin near Thessaloniki, Greece. The wider horizons opened by the described computational tools imply the need for detailed velocity structures, so to explain their impact on the radiated broad-band wavefield. Historically, long-period synthetic waveforms for large seismic events made use of the empirical Green’s function method, based on the self-similarity principle (Aki 1967) and on the representation theorem (Spudich & Archuleta 1987), and firstly proposed by Hartzell (1978) and Irikura (1983). Empirical approaches are however suitable just in the favourable yet rare condition of a vast number of recordings available for the site of interest or a very similar one (Chaljub et al.2015). However, a general concept was mediated from the mentioned traditional approach, that is, the idea that seismic observations recorded during small aftershocks identify with the system response to point-wise impulsive source excitation. From this point of view, the simulation of small aftershocks clarifies the impact of the propagating path and site effects, whose nature is manifold and sometimes counter-intuitive. One of the most investigated and debated aspects is the effect of the buried topography on the surface wavefield, due to the presence of complex geological interfaces such as basins, valleys, foldings and large heterogeneous soil deposits. The latter issues, along with the surface topography are responsible of wave scattering and diffraction, of body wave (de-)amplification, focusing and multiple reflection and sometimes of the presence of long-period surface waves. These effects may amplify due to proximity of the source to the site and, from a structural engineering point of view, they translate into a spatially incoherent incident motion, which entails the need for more refined seismic analysis when critical structures (such a nuclear power plants) are designed. For instance, Paolucci et al. (2015) recently provided an interesting insight on the role of the folded buried structure underneath the Po Plain, as one of the prominent causes of the ground motion distribution observed during the 2012 Emilia earthquake (Italy). Further interesting studies on the influence of the 3-D geological configurations on the characteristics of the ground shaking may be found in (Bard & Bouchon 1980; Faccioli & Vanini 2003; Gélis & Bonilla 2012; Chaljub et al.2015,among others). In this paper, the impact of a folded geology on the predicted ground shaking scenario of the site response of a nuclear power plant is assessed. The reference test case is represented by the 2007 MW6.6 Niigata-Ken Chūetsu-Oki earthquake (NCOEQ-2007), which struck a 100 km-wide area along the Japanese mid-west coastal area (depicted in Fig. 1), encompassing the Kashiwazaki-Kariwa Nuclear Power Plant (KKNPP). Figure 1. View largeDownload slide Map of the Niigata region, surrounding the Kashiwazaki-Kariwa Nuclear Power Plant (black square). The NCOEQ-2007 main shock is indicated as a red star, whereas coloured squares indicate the relocated aftershocks positions (TEPCO 2007). ASP1, ASP2 and ASP3 indicate the three major asperities deducted by Shiba (2008), by means of waveform inversion. Figure 1. View largeDownload slide Map of the Niigata region, surrounding the Kashiwazaki-Kariwa Nuclear Power Plant (black square). The NCOEQ-2007 main shock is indicated as a red star, whereas coloured squares indicate the relocated aftershocks positions (TEPCO 2007). ASP1, ASP2 and ASP3 indicate the three major asperities deducted by Shiba (2008), by means of waveform inversion. The KKNPP is the largest nuclear power plant in the world, with its seven units deployed around a 1.5 km wide area. The facility is located in the southern Niigata prefecture, adjacently the affected area (the Joyner-Boore and rupture distances were RJB = 0 km and RRUP = 16 km respectively; Yee et al.2011) and specifically on the hanging wall of the buried reverse fault which caused the earthquake. It has been reported that the site experienced nearly twice the value of Peak Ground Acceleration (PGA) assumed for the plant design (Pavlenko & Irikura 2012), but even despite the large ground settlement occurred, the redundant and over-designed structures incurred into minor damage (Kayen et al.2007). Nevertheless, in this paper, the spatial variability of the ground motion observed within the site area is investigated, based on previous studies which attributed it to the presence of a geological folding structure near the KKNPP (Tokumitsu et al.2009; Watanabe et al.2009; Hijikata et al.2010). For instance, Watanabe et al. (2009) and Watanabe et al. (2011) performed some 2-D FEM analyses, reliable up to 5.0 Hz, considering two cross-sectional of the mentioned folded structure. They investigated and successfully reproduced the amplification factors at two sites within the KKNPP, picturing the higher amplification occurred in the southernmost area. Following their work, Tsuda et al. (2011) extruded the 2-D cross sections along its directrix (approximately following the NW-SE coastline) and performed an FDM simulation of the NCOEQ-2007 main shock and three aftershocks, depicting the regional wavefield and the effect of an extended fault rupture, up to a frequency of 4.0 Hz. Aochi et al. (2013) tested three different 3-D geological structures (by simulating two aftershocks in the NCOEQ-2007 cluster) in an FDM simulation to reproduce the main shock scenario at long periods (i.e. up to 0.5 Hz). Aochi & Yoshimi (2016) improved the latter computational model, including the best available geological model and testing two different source models (both kinematic and dynamic descriptions, in the 0.1–1.0 Hz frequency band). Further pieces of research concerning the 3-D numerical simulations of the NCOEQ-2007 earthquake are the one provided by Kawabe & Kamae (2010) and Quinay et al. (2013). The formers simulated the wave propagation in a 3-D structure model, in a frequency range 0.05–1.6 Hz. Quinay et al. (2013) performed instead an entire fault-to-structure analysis, including the crustal geology, the soil strata and the KKNPP buildings. They tackled and solved this cumbersome computational problem by means of a two-step multiscale analysis based on (1) a solution of a low-resolution model (at a regional scale), used as an input boundary condition to the higher-resolution one (at an engineering-length scale). The computational scheme (formulated by means of the FEM) implied a domain decomposition with pre-partitioning, solved on parallel architecture. The authors obtained realistic results, compared to the observations in the Niigata region and at the KKNPP, being restrained to the frequency domain 0.1–1.0 Hz. In this context, most of the studies targeting the regional wavefield and the seismic source characterization lack of poor frequency resolution (see for instance Aochi & Yoshimi 2016; Quinay et al.2013), whereas models with higher frequency ranges (such as Watanabe et al.2009; Tsuda et al.2011) are focalized on the KKNPP subsurface structure. In this work, a 3-D regional-scale earthquake scenario of the area was built-up for two of the aftershocks belonging to the NCOEQ-2007 sequence and indicated as blue and green stars respectively in Fig. 1. Section 2 describes the main aspects of the observed waveforms at the KKNPP, focusing on the impulsive nature of the recorded velocigrams and the incoherent amplification between the southwest and northeast part of the site, placed at approximately 1 km from each other. The details on the geological structure in the Niigata region are provided in Section 3, stressing the methodological approach followed to conceive a fancy geological model suitable for reliable broad-band earthquake ground-motion simulations. The selection of the coarser and horizontally layered geology is described in Section 4 (the so called model LARGE), with the aim at reproducing the regional wavefield (solved by means of the Spectral Element Method, SEM) in a rather low-frequency range (up to fmax = 3.75 Hz, given the lack of a more detailed information at this scale) but constrained by several recording stations in the Niigata region. The impact of the folding is assessed in a second phase (Section 5) by plugging its 3-D configuration into the previous regional model (and obtaining the SEM model called SMALL). The outcome of numerical modelling from two aftershocks occurred at opposite sides of the KKNPP is analysed, comparing the synthetics (filtered at 5.0 Hz) at different locations within the nuclear plant. Special attention is given to the differences in terms of response spectra (Sa) at the ground surface, highlighting the discrepancies at different locations (corresponding to different Units of the nuclear site) and to their variation according to the aftershock position with respect to KKNPP. The comparison with the observations stress the need for a more accurate prediction of incident wave motion to be used as input for structural analyses, dropping the simplified assumption of layered half-space and supporting the importance of the complex geology as responsible of the ground motion incoherence at surface. 2 THE NCO EARTHQUAKE: OBSERVED STRONG GROUND MOTION CHARACTERISTICS The KKNPP site consists of seven Units, grouped into two blocks: (1) Units 1–4, built in the southwest part of the site and (2) Units 5–7 placed approximately 1.5 km away, in the northeast corner. Fig. 2 shows a sketch of the plant, highlighting the recording devices downhole and at surface. Devices 1G1 (green) and 5G1 (orange) both belong to the most recent (2004) recording network, and they are placed at the surface (G.L. 0 m), respectively close by Unit 1 and Unit 5–7. On the other hand, four devices (SG1–SG4) are deployed along the vertical array KSH (blue) located at the KKNPP Service Hall, up to a depth of 250 m; other five devices are deployed downhole close by Unit 5 (KK5, red), till 312 m of depth; four devices are finally located in the surroundings of Unit 1 (KK1, magenta) with devices G07–G10 reaching 255 m of depth. Figure 2. View largeDownload slide Map of the KKNPP site (courtesy of TEPCO 2007). The five coloured squares indicate the three recording stations downhole (KSH at Service Hall (blue), KK1 at Unit 1 (magenta), KK5 at Unit 5 (red)) and the two surface ones free-field (1G1 at Unit 1 (green) and 5G1 at Unit 5 (orange)). The KSH station downhole, along with recording devices at surface 1G1 and 5G1 entirely recorded the NCOEQ-2007 main shock. The devices were oriented with respect to plant North (NS-KKNPP), which differs from the real geodetic North of an angle θ = 18° 54$$^{^{\prime }}$$ 51$$^{^{\prime \prime }}$$. Moreover, some of them suffered of an azimuthal deviation, whose value was provided by TEPCO (2007). Figure 2. View largeDownload slide Map of the KKNPP site (courtesy of TEPCO 2007). The five coloured squares indicate the three recording stations downhole (KSH at Service Hall (blue), KK1 at Unit 1 (magenta), KK5 at Unit 5 (red)) and the two surface ones free-field (1G1 at Unit 1 (green) and 5G1 at Unit 5 (orange)). The KSH station downhole, along with recording devices at surface 1G1 and 5G1 entirely recorded the NCOEQ-2007 main shock. The devices were oriented with respect to plant North (NS-KKNPP), which differs from the real geodetic North of an angle θ = 18° 54$$^{^{\prime }}$$ 51$$^{^{\prime \prime }}$$. Moreover, some of them suffered of an azimuthal deviation, whose value was provided by TEPCO (2007). The three vertical arrays KSH, KK5 and KK1 belong to the oldest recording network, where surface devices 1G1 (Unit 1, G.L. 0 m) and 5G1 (Unit 5, G.L. 0 m) are deployed at surface. It is worth noting that just one (KSH, Service Hall) out of three vertical array entirely recorded the main shock, along with 1G1 and 5G1. According to Watanabe et al. (2009), the strong motion records at Unit 1 were significantly larger than those in the surrounding of Unit 5. Many authors (e.g. Uetake et al.2008; Miyake et al.2010; Tsuda et al.2011) stressed the impulsive nature of the recorded time-histories. Large velocity pulses are usually observed nearby the fault rupture and coupled with their corresponding large peak displacement, they considerably enhance the structural damage potential (Cox & Ashford 2002). In the NCOEQ-2007 case, Uetake et al. (2008) first noticed the three significant pulses at KSH array, associating them to the three major asperities identified (via waveform inversion) on the fault plane. As a matter of fact, velocity pulses may be effectively extracted from the time-histories at KSH (SG1, G.L. −2.4 m and SG4, G.L. −250 m), 1G1 and 5G1 (i.e. the only devices which entirely recorded the main shock), by employing the ranking criterion proposed by Baker (2007) (excluding late arrivals and small events). Fig. 3 shows an example of the mentioned pulse-like nature of the velocigrams recorded. Figure 3. View largeDownload slide Comparison between velocity signal (black traces) and extracted pulses (red traces), according to Baker (2007), for devices (a) SG1 (G.L. −2.4 m), (b) SG4 (G.L. −250 m) in the KSH vertical array (Service Hall) and for (c) 1G1 (G.L. 0 m) and (d) 5G1 (G.L. 0 m). KKZ indicates the set of new seismometers, installed in 2004, to which 1G1 and 5G1 belong. The in plane azimuthal angle θ at which those pulses were extracted is indicated above the velocigrams. Figure 3. View largeDownload slide Comparison between velocity signal (black traces) and extracted pulses (red traces), according to Baker (2007), for devices (a) SG1 (G.L. −2.4 m), (b) SG4 (G.L. −250 m) in the KSH vertical array (Service Hall) and for (c) 1G1 (G.L. 0 m) and (d) 5G1 (G.L. 0 m). KKZ indicates the set of new seismometers, installed in 2004, to which 1G1 and 5G1 belong. The in plane azimuthal angle θ at which those pulses were extracted is indicated above the velocigrams. The average pulse period ranges around 2.5 s for SG1, SG4 and 5G1, whereas 1G1 identifies with a shorter period pulse (1.5 s). Remarkable differences (in terms of amplitude, period and phase) are evident at a first glance, seemingly indicating an incoherency of the ground motion all around the nuclear facility zone. The observed site response at the surface is definitely more intense at the Service Hall and at Unit 1. It is worth noting that the velocity pulses were identified and extracted at different azimuthal directions (with respect to the real North) for SG1 (KSH, surface) SG4 (KSH, depth), 1G1 (Unit 1, surface) and 5G1 (Unit 5, surface) respectively, the polarization of the impulsive incident wavefield spatially varies at the site scale, along its travelling path towards the surface. Another interesting aspect resides in the fact that later pulses, whose nature has been correlated to the progressive fault rupture of the three asperities, may be identified at each location (Uetake et al.2008; Miyake et al.2010). Tsuda et al. (2011) observed that the first two pulses were rather comparable in terms of maximum amplitude, whereas the third one (allegedly coming from the third asperity ASP3, represented as green rectangle in Fig. 1) is different for two sites: the peak ground velocity on the south side of the KKNPP (KK1, downhole array at Unit 1) is much larger than that on the north side (KK5, downhole array at Unit 5) (Tsuda et al.2011). 3 CONSTRUCTION OF THE HYBRID GEOLOGICAL MODEL The subsurface geological structure underlying the Niigata region has been proven to be rather intricate. Several 1-D and 3-D models of the velocity structure have been proposed in the past, based on geological and geophysical explorations (see, for instance, Shinohara et al.2008). The latter set of geological configurations have been tested and tuned by several authors, by exploiting the long aftershock sequence occurred in the main shock aftermath and then exploited to simulate the main shock itself. Table 1 resumes the ensemble of ground shaking analyses performed on the KKNPP site, during the NCOEQ-2007. Kawabe & Kamae (2010) simulated the wave propagation in a 3-D structure model provided by the Japan Nuclear Energy Safety Organization (JNES, internal report, 2005) in a frequency range 0.05–1.6 Hz. Aochi et al. (2013) and Aochi & Yoshimi (2016) performed a complete FDM analysis of the NCOEQ-2007 scenario, comparing three regional velocity models: (1) the model proposed by the Earthquake Research Institute (ERI; University of Tokyo) obtained by P- and S-wave traveltime double-difference tomography on the data provided by Shinohara et al. (2008); Kato et al. (2008, 2009), with a grid resolution of 3 km × 5 km × 3 km (vertical) and with a minimum shear wave velocity VS of 866 m s−1; (2) the model proposed by the National Research Institute for Earth Science and Disaster Resilience (NIED), taken from Fujiwara et al. (2009) from the seismic reflection results and observed H/V spectra (available on the Japan Seismic Hazard Information Station (J-SHIS), website: http://www.j-shis.bosai.go.jp), with a 1 km × 1 km resolution and minimum VS equal to 350 m s−1; (3) the model proposed by the Geophysical Survey of Japan (GSJ) and found in Sekiguchi et al. (2009), which is an improvement of the NIED velocity structure for the Niigata area (the resolution is 0.5 km × 0.5 km), carefully tuned upon the regional variation of material parameters (VSmin = 400 m s−1). The authors concluded that the ERI model was less preferable than the other two in order to reproduce the regional wave field, mainly due to the poor resolution of the shallow deposits and despite the low frequency range of the analyses. Moreover, they argue the use of the JNES model (Kawabe & Kamae 2010) for the regional-wise purposes, being however more precise around KKNPP. The crust material, topography, and subsurface data provided by the J-SHIS were used by Quinay et al. (2013) to carry out a complete 3-D analysis up to 1.0 Hz, studying the effect of the synthetic wave field on the structural behaviour of one of KKNPP reactor building. Aochi & Yoshimi (2016) improved the previous results by Aochi et al. (2013) by enlarging the frequency band up to 1.0 Hz and employing the GSJ model, although their study targeted mostly the effect of different seismic source models and numerical implementation. Table 1. Summary of ground shaking earthquake ground shaking scenarios of the NCOEQ-2007, targeting the KKNPP. FDM, Finite Difference Method; FEM, Finite Element Method; SEM, Spectral Element Method. Fold, Tokumitsu et al. (2009); JNES-2005, internal report; ERI, Shinohara et al. (2008); Kato et al. (2008, 2009); JNES-2008, JNES (2008); NIED, Fujiwara et al. (2009); GSJ, Sekiguchi et al. (2009); Reference  Dimension  Geology  Method  Model Size [km]  fmax [Hz]  min(VS) $$[\rm {{m}\,{s}}^{-1}]$$  Watanabe et al. (2009)  2-D  Fold  FEM  7.6 × 4.8  5.0  700  Kawabe & Kamae (2010)  3-D  JNES-2005  ?  ?  1.6  700  Ducellier & Aochi (2010)  3-D  ERI  FDM  110 × 120 × 30  0.866  866  Tsuda et al. (2011)  3-D  JNES-2008 + Fold  FDM  50 × 50 × 20  4.0  700  Aochi et al. (2013)  3-D  ERI/NIED/GSJ  FDM  110 × 120 × 30  0.5  866/350/400  Quinay et al. (2013)  3-D  NIED  FEM  110 × 120 × 30  1.0  866/350/400  Aochi & Yoshimi (2016)  3-D  GSJ  FDM  110 × 120 × 30  1.0  400  LARGE (this work)  3-D  Aochi2013  SEM  90 × 83 × 82  3.75  1020  SMALL (this work)  3-D  Aochi2013+Fold  SEM  68 × 50.8 × 50.8  5.0  700  Reference  Dimension  Geology  Method  Model Size [km]  fmax [Hz]  min(VS) $$[\rm {{m}\,{s}}^{-1}]$$  Watanabe et al. (2009)  2-D  Fold  FEM  7.6 × 4.8  5.0  700  Kawabe & Kamae (2010)  3-D  JNES-2005  ?  ?  1.6  700  Ducellier & Aochi (2010)  3-D  ERI  FDM  110 × 120 × 30  0.866  866  Tsuda et al. (2011)  3-D  JNES-2008 + Fold  FDM  50 × 50 × 20  4.0  700  Aochi et al. (2013)  3-D  ERI/NIED/GSJ  FDM  110 × 120 × 30  0.5  866/350/400  Quinay et al. (2013)  3-D  NIED  FEM  110 × 120 × 30  1.0  866/350/400  Aochi & Yoshimi (2016)  3-D  GSJ  FDM  110 × 120 × 30  1.0  400  LARGE (this work)  3-D  Aochi2013  SEM  90 × 83 × 82  3.75  1020  SMALL (this work)  3-D  Aochi2013+Fold  SEM  68 × 50.8 × 50.8  5.0  700  View Large 3.1 Description of the layered geology of the Niigata region As a general remark, the lowest shear-wave velocity varies from 350 m s−1 to 866 m s−1, for the studies previously, estimated for low-frequency analyses (i.e. ground motion reproduction and a few waveform inversions of the source mechanism), and eventually based on the interpretation of several 1-D/2-D models, which Aochi et al. (2013) proved to perform fairly well in reproducing the complex wavefield radiated, although phase shifts must be taken into account, due to modified traveltimes. Fig. 4 shows some of the mentioned layered velocity structures available for the Niigata region. D&A2010-1 and D&A2010-2 were proposed by Ducellier & Aochi (2010), in their numerical simulations performed for the NCOEQ-2007 main shock. With the same purposes, Aochi et al. (2013) provided a simplified 1-D profile, tagged as Aochi2013. Cire.2008 velocity values were employed by Cirella et al. (2008) to perform a waveform inversion of the NCOEQ-2007 source mechanism, whereas Shino.2008 was obtained by exploiting the data of the seismic survey conducted by Shinohara et al. (2008) and used to relocate the aftershock sequence. Finally, profile indicated as NIG004, NIG016, NIG026 and NIGH12 refer to the namesake K-NET and KiK-net stations. Figure 4. View largeDownload slide Set of regional 1-D velocity structures (VP and VS) for the Niigata region, estimated by several authors. D&A2010-1 and D&A2010-2 were proposed by Ducellier & Aochi (2010). Aochi et al. (2013) and Cirella et al. (2008) provided the profiles tagged as Aochi2013 and Cire.2008 respectively. Shino.2008 was obtained by Shinohara et al. (2008) and profiles indicated as NIG004, NIG016, NIG026 and NIGH12 refer to the namesake K-NET and KiK-net stations (the NIED strong ground motion networks, covering the Japan territory) and they were estimated by K. Koketsu (source: http://taro.eri.u-tokyo.ac.jp/saigai/chuetsuoki/source/index.html) Figure 4. View largeDownload slide Set of regional 1-D velocity structures (VP and VS) for the Niigata region, estimated by several authors. D&A2010-1 and D&A2010-2 were proposed by Ducellier & Aochi (2010). Aochi et al. (2013) and Cirella et al. (2008) provided the profiles tagged as Aochi2013 and Cire.2008 respectively. Shino.2008 was obtained by Shinohara et al. (2008) and profiles indicated as NIG004, NIG016, NIG026 and NIGH12 refer to the namesake K-NET and KiK-net stations (the NIED strong ground motion networks, covering the Japan territory) and they were estimated by K. Koketsu (source: http://taro.eri.u-tokyo.ac.jp/saigai/chuetsuoki/source/index.html) 3.2 Description of the local 3-D geology Despite the good approximation of the regional incident wave motion obtained with simplified 1-D and 3-D geological configurations in the surroundings of the epicentral area, there is consensus among many researchers (Uetake et al.2008; Tokumitsu et al.2009; Watanabe et al.2009; Hijikata et al.2010) that none of them is suitable to picture the spatial variability of the ground motion observed at KKNPP. With this regard, Uetake et al. (2008) and Tokumitsu et al. (2009) related this aspect to the presence of a folded structure underneath KKNPP (characterized via boring and seismic reflection surveys, Kobayashi et al.1995) and composed by the Madonosaka syncline, interposed between the Ushirodani and Chuo-Yutai anticlines (Gürpinar et al.2017). The folding cross-section extends for approximately 7 km across the Japanese coastline, up to 2.5 km down depth and with its hinge axes (i.e. the directrixes) strike at N55°E, as shown in Figs 5(a) and (b). The geographical distribution of this peculiar sediment conformation shears through the nuclear facility, with the Ushirodani anticline and Madonosaka syncline placed below Unit 5 and Unit 1 respectively. The Chuo-Yatai anticline is indicated as well, which seemingly passes in the Service Hall surroundings (array KSH). A detailed description of the seven strata defining the folding is provided in Section 5. Watanabe et al. (2009) constructed a 2-D FEM model of the syncline-anticline conformation, from boring and seismic reflection survey. Simulation results show good agreement with the observed strong motion records at Unit 1. Tsuda et al. (2011) constructed a 3-D model by taking into consideration seven 2-D sections across the folding area, and interpolating them. The joined up the folding model with the regional one (the 3-D velocity model built by the Japan Nuclear Energy Safety Organization (JNES), which covers broad area including Chuetsu area (JNES 2008). 3.3 Proposed hybrid model of the KKNPP earthquake scenario Given the complex geological structure described so far, and with the twofold aim to reproduce the overall regional wavefield, yet focusing on a realistic broad-band simulation of the incident ground motion in the KKNPP surroundings, a hybrid geological model of the Niigata region was herein constructed in two steps, progressively adding all the available information. First of all, the layered Aochi2013 profile (described in Section 3.1) was selected as representative for the propagation of the long period ground motion (i.e. 0.1–3.75 Hz) at a regional scale (i.e. approximately 90 km × 83 km × 82 km). From now on, it will be referred as to the LARGE model. Its considerable extension was mainly chosen to grant a minimal coverage of the NW regional wavefield. Thus LARGE model was designed so as to include part of the Sado Island in front of the Niigata region, where the K-NET station NIG004 is located. The Japan sea presence was disregarded at this step, along with related bathymetry and coastline (see Fig. 5b). Figure 5. View largeDownload slide (a) Schematic structural map of the series of Ushirodani anticline—Madonosaka syncline, located underneath KKNPP. Red, orange and fuchsia points indicate the KKNPP site and AS1 and AS2 epicentres, respectively. The white line represents section SC. (b) Sketch of the extension of the computational models (LARGE and SMALL) considered in this study. Figure 5. View largeDownload slide (a) Schematic structural map of the series of Ushirodani anticline—Madonosaka syncline, located underneath KKNPP. Red, orange and fuchsia points indicate the KKNPP site and AS1 and AS2 epicentres, respectively. The white line represents section SC. (b) Sketch of the extension of the computational models (LARGE and SMALL) considered in this study. Section 4 presents the calibration of LARGE model, dealing with the choice of an appropriate source time function and topographical effect. The LARGE model was trimmed and adjusted so as to plug the 3-D Ushirodani anticline—Madonosaka syncline—Chuo-Yatai anticline structure (widely described in Section 3.2) into it, following the previous work of Tsuda et al. (2011) (see Section 5). The 1-D profile, issued from Aochi2013, was welded to the syncline-anticline system by adjusting and linearly smoothing the uppermost crustal material discontinuities to reconnect with the external layered geology of LARGE model. 4 SIMULATION OF REGIONAL WAVEFIELD: LARGE MODEL 4.1 Preliminary selection of a suitable geological profile The epistemic uncertainty on the regional geology was outlined and solved by first employing the Wave-Number-Integration method (WNI), proposed and implemented by Hisada (1994, 1995, 2008), that is, a semi-analytical method that simulates the complete 3-D wave propagation field (based on the computation of static and dynamic Green’s functions of displacements and stresses) radiated from an extended kinematic seismic source in a viscoelastic horizontally layered half-space. The MJMA4.4 aftershock of July 16, 21.08h (JST), was considered, which nucleated near-by one of the three major fault asperities at a depth of 11 km (Tsuda et al.2011) determined by detecting the polarity for the P-wave first motion (performed by the Hi-net (Obara et al.2005)) and by considering the Centroid Moment Tensor analysis provided by F-net data (Kubo et al.2002). The KiK-net and K-NET recording stations were employed at this point, along with the recorded time-histories at the KKNPP (device KSH-SG4, placed at G.L.-250m at the Service Hall of the nuclear site, see Fig. 2). Several WNI analyses were performed, testing the estimated geology in Fig. 4. This preliminary investigation ended up with the choice of the Aochi2013 soil profile proposed by Aochi et al. (2013) and whose velocity and attenuation properties are summarized in Table 2. Fig. 6 shows a fairly good comparison obtained at several stations in the surroundings of the KKNPP. In accordance with Aochi et al. (2013), the recordings and the synthetics were band-pass filtered between 0.1–0.5 Hz and aligned at the time-step corresponding to the 1 per cent of the respective Arias intensity. Although a satisfactory match has been obtained for other soil profiles (in fact, similar results were obtained for the crustal models tagged as Cire.2008 and Shin.2008 in Fig. 4), Aochi et al. (2013) introduced a sequence of softer layers (Fig. 4) at shallow depths, with respect to the crustal models estimated for the K-NET/KiK-net stations (NIG004, NIG016, NIG026, NIGH12). Despite the local alignment, some phase shift among the simulations and the recordings is still visible for NIG004 and NIGH12, whereas in the near field (cf. NIG018, NIG016) the fit is quite good (probably due to the simple propagation path, Aochi et al.2013)). Finally, from Fig. 6, it appears that the simulations performed at KKNPP site are satisfactory only along the NS direction, whereas synthetics are de-amplified along the EW direction. This might be an effect of the shallow local geology of the Niigata basin. With that being said, it is well known that the choice of a 1-D velocity structure tends to accentuate the ground motion directionality. Figure 6. View largeDownload slide Simulations performed by WNI (red line) method of the NCOEQ-2007 MJMA4.4 aftershock of July 16, 21.08h (JST). Blue waveforms (REC) represent the recorded velocigrams, red waveforms (WNI) the synthetic ones. Both records and synthetics were base-line corrected and bandpass filtered between 0.1 and 0.5 Hz. Synthetic waveforms were obtained by considering the soil profile Aochi2013 (Fig. 4). Velocigrams (in cm s−1) are herein magnified by a factor 1000 and aligned with respect to the 1 per cent of their Arias intensity. Figure 6. View largeDownload slide Simulations performed by WNI (red line) method of the NCOEQ-2007 MJMA4.4 aftershock of July 16, 21.08h (JST). Blue waveforms (REC) represent the recorded velocigrams, red waveforms (WNI) the synthetic ones. Both records and synthetics were base-line corrected and bandpass filtered between 0.1 and 0.5 Hz. Synthetic waveforms were obtained by considering the soil profile Aochi2013 (Fig. 4). Velocigrams (in cm s−1) are herein magnified by a factor 1000 and aligned with respect to the 1 per cent of their Arias intensity. Table 2. Properties of the Aochi2013 profile (Aochi et al.2013). z represents the depth of the upper layer surface, VP and VS are the P-wave and S-wave velocities respectively, QP and QS the quality factors for P-wave and S-wave respectively. The ⋆ indicates the interface chosen to plug the folding structure into the 1-D Aochi2013 profile, granting a smooth transition from one model to the other. z[km]  $$V_P [\rm {{km}\,{s}}^{-1}]$$  QP[1]  $${V}_{S} [\rm {{km}\,{s}}^{-1}]$$  QS[1]  0.0  2.28  200.0  1.02  100.0  0.5  2.57  228.7  1.23  100.0  1.0  2.93  257.1  1.51  100.0  1.5  3.09  293.3  1.63  100.0  2.0  3.25  309.1  2.09  100.0  3.0  3.68  325.9  2.33  100.0  4.0  4.03  368.2  2.49  100.0  ⋆5.0  4.30  403.2  2.75  100.0  6.0  4.55  430.9  2.89  100.0  7.0  4.76  455.8  3.10  100.0  8.0  5.00  476.4  3.40  100.0  9.0  5.37  500.2  3.76  100.0  10.0  5.88  537.5  3.80  100.0  14.0  6.51  588.7  3.85  250.0  z[km]  $$V_P [\rm {{km}\,{s}}^{-1}]$$  QP[1]  $${V}_{S} [\rm {{km}\,{s}}^{-1}]$$  QS[1]  0.0  2.28  200.0  1.02  100.0  0.5  2.57  228.7  1.23  100.0  1.0  2.93  257.1  1.51  100.0  1.5  3.09  293.3  1.63  100.0  2.0  3.25  309.1  2.09  100.0  3.0  3.68  325.9  2.33  100.0  4.0  4.03  368.2  2.49  100.0  ⋆5.0  4.30  403.2  2.75  100.0  6.0  4.55  430.9  2.89  100.0  7.0  4.76  455.8  3.10  100.0  8.0  5.00  476.4  3.40  100.0  9.0  5.37  500.2  3.76  100.0  10.0  5.88  537.5  3.80  100.0  14.0  6.51  588.7  3.85  250.0  View Large 4.2 Verification of the model upgrade (WNI→SEM) A Spectral Element Method software, SEM3D, was used to solve the 3-D wave propagation problem so as to verify the code with the WNI semi-analytical model tuned in the earlier stages of the parametric analysis. The chosen 1-D geological model was implemented in 3-D large scale SEM numerical model. No topography was introduced at first, so as to compare the synthetics with WNI results. Five Gauss–Lobatto–Legendre (GLL) integration nodes (the integration points for the SEM formulation) were chosen to grant the minimum wavelength (≈500 m) a satisfactory discretization (i.e. to reach a frequency upper bound of ≈2 Hz, with a third order polynomial degree). In the SEM formulation, material properties are assigned to each GLL point. Therefore, two strategies were investigated to feature the numerical model with the varying material properties, namely (1) a direct meshing approach that adapts the computational grid to the layer-to-layer interface, so to coherently reproduce the abrupt impedance contrast along it and (2) a not honouring approach that associate to the whole independently meshed domain a space-varying heterogeneous material map (in this case, transverse isotropic in the horizontal plane). In this standard test case (i.e. a flat and bounded horizontally layered half-space) the second strategy does not seem particularly appealing since the material interfaces are planar and horizontal and a structured mesh was adopted. The numerical verification of the model upgrade (WNI→SEM3D) was satisfactory. For instance, Fig. 7 shows a synoptic comparison between synthetic waveforms filtered between 0.1 and 1.0 Hz. Although no computational gain is achieved by using one or the other approach at this early stage, the verification test was performed in view of implementing the not honouring approach in the non-structured mesh domain, extruded at depth from the local Digital Elevation Model. Figure 7. View largeDownload slide Numerical verification of SEM3D code, against the semi-analytical solution obtained via WNI method (0.1-1.0 Hz). Synthetic velocity waveforms obtained at (a) the KKNPP site (KSH-SG4, G.L. −250 m, focal distance 15.1 km) and at (b) the K-NET station NIG018 (focal distance of 20.1 km) are portrayed. A good match is achieved between the WNI method (WNI, red line) and SEM3D numerical model equipped either with a layered geology (LAY, blue line) or with a not honouring approach (GRD, green line) Figure 7. View largeDownload slide Numerical verification of SEM3D code, against the semi-analytical solution obtained via WNI method (0.1-1.0 Hz). Synthetic velocity waveforms obtained at (a) the KKNPP site (KSH-SG4, G.L. −250 m, focal distance 15.1 km) and at (b) the K-NET station NIG018 (focal distance of 20.1 km) are portrayed. A good match is achieved between the WNI method (WNI, red line) and SEM3D numerical model equipped either with a layered geology (LAY, blue line) or with a not honouring approach (GRD, green line) 4.3 Effect of the surface topography A further step of the numerical exercise consisted into adding the surface topography (TOPO) to the previously calibrated flat model (FLAT). 15 Gauss-Lobatto-Legendre (GLL) integration nodes were chosen to propagate a seismic wavefield up to ≈3.75 Hz (considering, as reference, a third-order polynomial interpolation, Paolucci et al.2015). The shear-wave velocity profile (Aochi2013 model) was implemented via a not honouring approach. The Sea of Japan is not taken into account: a flat solid surface was considered instead. Fig. 8 shows the comparison between the accelerograms obtained at the ground level at KKNPP site, for the FLAT (red line) and TOPO (blue line) models. The waveforms were band-pass filtered between 0.05 and 3.75 Hz. Early wave arrivals at higher frequencies are observed when the surface topography is included in the model, although it does not seem to consistently affect the radiated wavefield. This aspect might be due to the fact that the Japan sea was disregarded. Figure 8. View largeDownload slide Numerical verification of SEM3D code: effect of the surface topography. Synthetic waveforms obtained at KKNPP (at surface) are portrayed: (a) EW component, (b) NS component. Both synthetics were compared in the frequency band 0.05–3.75 Hz. Figure 8. View largeDownload slide Numerical verification of SEM3D code: effect of the surface topography. Synthetic waveforms obtained at KKNPP (at surface) are portrayed: (a) EW component, (b) NS component. Both synthetics were compared in the frequency band 0.05–3.75 Hz. 4.4 Comparison with ground motion prediction equations The a priori selection of a suitable Source Time Function (STF) for a point-wise double-couple approximation of the seismic source is not obvious (Bizzarri 2014). To clarify this fundamental aspect, a brief parametric analysis addressing different STFs and the respective featuring parameters was performed. For instance, Aochi et al. (2013) simulated two aftershocks as double-couple point-wise sources, with smooth bell-shaped STF (i.e. cubic B-spline function) of a 0.5 s duration at the hypocentre. Specifically, two classical STF are addressed, namely (1) the so called Bouchon’s ramp (or smoothed ramp function) $$u_S^{\text{TNH}}\left(t\right)$$ and (2) the source model proposed by Brune (1970) $$u_S^{\text{EXP}}\left(t\right)$$. The former is well suited to approximate long-tailed Moment Rate Functions (MRF; Duputel et al.2013). Two are the free parameters that tune the shape and the position in time of the two mentioned STFs: the maximum displacement offset umax and the generic rise-time τR. The maximum slip value umax is spontaneously determined by relating it to the scalar seismic moment value (for the MJMA4.4 aftershock event, the F-net indicated a seismic scalar moment M0 = 5.21× 1015 Nm), whereas τR is rather delicate to be estimated. Although no explicit definition of the rise-time is available neither for $$u_S^{\text{TNH}}\left(t\right)$$ nor for $$u_S^{\text{EXP}}\left(t\right)$$ (Bizzarri 2014), τR represents the duration of the slip rate pulse (Bizzarri 2014; De Martin et al.2007). Whilst estimated fairly easily for strong ground shaking, empirical correlations may provide longer rise-times, due to the fact that they are based upon registrations of larger earthquake events, polluted by the effect of the rupture time along the fault plane. To assess the influence of τR on the radiated regional-wave field, three sources were tested on the LARGE model (reported in Table 3). Two values of τR were chosen, being τR1 the value of the rise-time proposed in Dreger et al. (2007), estimated as $$\tau _{R1} = 10^{0.5\left({\rm M}_{\rm W}-6.69\right)}$$, and τR2 being derived from the expression of the semi-duration proposed by Duputel et al. (2012) $$\tau _{R2}/2 = 1.2\cdot 10^{-8} \root 3 \of {{\rm M}_0\left[{\rm dyne.cm}\right]}$$. The latter value is closer to the rise-time indicated by Tsuda et al. (2011), based on the width of S-wave pulse of the observed waveforms at rock sites of F-net stations adjacent to the source (and equal to 0.7 s). Table 3. Source Time Functions (STFs) chosen to represent the point-wise seismic moment time-history. Brune’s pulse $$u_S^{{\rm EXP}}$$ refers to the STF proposed by Brune (1970), whereas the Bouchon’s ramp $$u_S^{{\rm TNH}}$$ refers to the STF proposed by Bouchon (1981). H(t) represent the Heaviside function. STF  Function  $$u_S^{{\rm EXP}}$$  $$u_{{\rm max} \left[1-\left(1+\frac{4t}{\tau _R}\right)\exp \left(-\frac{4t}{\tau _R}\right)\right]H\left(t\right)}$$  $$u_S^{{\rm TNH}}$$  $$\frac{u_{\rm max}}{2} \left[1+\tanh \left(\frac{4t}{\tau _R}\right)\right]$$  STF  Function  $$u_S^{{\rm EXP}}$$  $$u_{{\rm max} \left[1-\left(1+\frac{4t}{\tau _R}\right)\exp \left(-\frac{4t}{\tau _R}\right)\right]H\left(t\right)}$$  $$u_S^{{\rm TNH}}$$  $$\frac{u_{\rm max}}{2} \left[1+\tanh \left(\frac{4t}{\tau _R}\right)\right]$$  View Large The outcome of the parametric analysis is presented herein in terms of comparison with a classical Ground Motion Prediction Equation (GMPE), proposed by Fukushima (2007). The choice of this GMPE is justified by the coherency between the predicted values and aftershock recordings in the Niigata region, although for small magnitude events GMPEs are less controlled by the rupture process, being related more to stress drop and wave propagation effects. Fig. 9 portrays the mentioned comparison in terms of geometric mean of the horizontal PGAs and PGVs. Specifically, Figs 9(a)–(c) pertain the comparison between recordings and the Fukushima GMPEs for PGA and PGV respectively, whereas Figs 9(b)–(d) refer to the SEM3D synthetics. Figure 9. View largeDownload slide Recorded (a) and synthetic (b) horizontal geometric mean of Peak Ground Acceleration (PGAH) obtained by SEM3D simulations of the MJMA4.4 aftershock event, compared to the GMPE proposed by Fukushima (2007). Recorded (c) and synthetic (d) horizontal geometric mean of Peak Ground Velocity (PGVH) obtained by SEM3D simulations of the MJMA4.4 aftershock event, compared to the GMPE proposed by Fukushima (2007). Red circles and blue triangles refer to Brune (uEXP) and a Bouchon’s (uTNH) functions respectively, with free-parameter τR1 = 0.113s; green triangles refer to a Bouchon’s ramp function with free-parameter τR2 = 0.820s instead. Thick black line represents the median PGAH/PGVH stemming from the chosen GMPE; dashed black lines (±σ) represent the standard deviations. Figure 9. View largeDownload slide Recorded (a) and synthetic (b) horizontal geometric mean of Peak Ground Acceleration (PGAH) obtained by SEM3D simulations of the MJMA4.4 aftershock event, compared to the GMPE proposed by Fukushima (2007). Recorded (c) and synthetic (d) horizontal geometric mean of Peak Ground Velocity (PGVH) obtained by SEM3D simulations of the MJMA4.4 aftershock event, compared to the GMPE proposed by Fukushima (2007). Red circles and blue triangles refer to Brune (uEXP) and a Bouchon’s (uTNH) functions respectively, with free-parameter τR1 = 0.113s; green triangles refer to a Bouchon’s ramp function with free-parameter τR2 = 0.820s instead. Thick black line represents the median PGAH/PGVH stemming from the chosen GMPE; dashed black lines (±σ) represent the standard deviations. The Fukushima (2007) GMPE predicts rather well both the recorded mean PGAs and PGVs (±σ), although a slight overestimation may be noticed at larger source-to-site distances. This discrepancy might be due to the magnitude (greater or equal to 5.0) to calibrate the GMPE parameters (Cotton et al.2008). However, Figs 9(b)–(d) suggest that the choice of a rise-time τR1 (for either $$u_S^{{\rm EXP}}$$ or $$u_S^{{\rm TNH}}$$) is the most consistent with the selected GMPE (Fukushima 2007), so it has been chosen hereafter as the suitable STF for the aftershock analyses. It is worth noticing that the choice of a constant Q factor might affect the peak values, diminishing their values compared to the attenuation relationship. 5 SIMULATION OF NEAR-SOURCE WAVEFIELD: SMALL MODEL At this point, the skeleton geology (LARGE model) was upgraded by obtaining an hybrid one (the 3-D structure referred hereafter as SMALL) so to insert the folded structure of the Earth’s crust underneath KKNPP (see the geographical map in Fig. 10a). The latter 3-D geology is defined by seven strata: a representative cross-section SC (see Fig. 10d, corresponding to section S4 in Tsuda et al.2011), passing through KKNPP, was chosen and extruded up to SW (south-westward) and SE (north-eastward). The folding shape is linearly smoothed (the so called smoothing bands in Fig. 10a) so to reconnect with the external subhorizontally layered geology defined in the previous section. This choice was justified by the possible spurious noise generated by the boundaries of the folding zone and polluting the synthetics. Figure 10. View largeDownload slide (a) Geographical map of the KKNPP surroundings, depicting folded cross-section SC, the smoothing bands and the extrusion edged, cross-sections SE and SW. (b and c) Sketched of the smoothed folded structure plugged into the SMALL model: Section SW (b) and Section SC (c) The colours indicate different layers. Figure 10. View largeDownload slide (a) Geographical map of the KKNPP surroundings, depicting folded cross-section SC, the smoothing bands and the extrusion edged, cross-sections SE and SW. (b and c) Sketched of the smoothed folded structure plugged into the SMALL model: Section SW (b) and Section SC (c) The colours indicate different layers. Table 4 summarizes the mechanical properties assumed in this paper for the folding cross-section SC (firstly proposed by Watanabe et al. (2009, 2011); Tsuda et al. (2011)). Table 4. Geological properties of the folding structure underneath KKNPP. VP and VS are the pressure- and shear-wave velocities respectively. The ⋆⋆ indicates the interface chosen to plug the folding structure into the original 1-D Aochi2013 profile, granting a smooth transition from one model to the other. Layer  VS[m s − 1]  VP[m s − 1]  ρ[Kg m − 3]  Nishiyama  700  1900  1700  Shiiya  1200  2200  2100  Upper Teradomari  1700  3300  2300  Lower Teradomari  2000  4200  2400  Nanatani  2000  4600  2500  ⋆⋆Green tuff  2600  5200  2600  Seismic bed rock  2600->2750  5200  2600  Layer  VS[m s − 1]  VP[m s − 1]  ρ[Kg m − 3]  Nishiyama  700  1900  1700  Shiiya  1200  2200  2100  Upper Teradomari  1700  3300  2300  Lower Teradomari  2000  4200  2400  Nanatani  2000  4600  2500  ⋆⋆Green tuff  2600  5200  2600  Seismic bed rock  2600->2750  5200  2600  View Large Although the mechanical properties of the folding model do not differ much from the Aochi2013 profile below 2 km deep, the former geological profile is featured by thinner strata and by lower minimum shear-wave velocity. However, a slight correction in the values of the seismic bedrock shear wave velocity VS was performed herein. The latter was adjusted to grant a smooth transition from the Aochi2013 horizontally layered regional structure to the folding detail. To naturally plug the latter into the former, two welding surfaces were selected, being respectively the upper surface of Aochi2013 layer located at 5 km depth (whose estimated VS value is 2750.0 m s−1, indicated with ⋆ in Table 2) and the lower interface of the Green tuff layer (whose estimated VS value is 2600.0 m s−1, indicated with ⋆⋆ in Table 4) as counterpart in the folding structure. In between the two surfaces, a linear variation of the VS, VP and ρ values is considered. Just replacing the local folding model within the broad Aochi2013 model would produce (as shown by Tsuda et al. (2011) artificial refracted waves, generated around the boundaries/gaps into the hybrid velocity structure. In the SMALL model, moreover, the Japan sea was considered. Due to the asymmetric nature of the folded structural geology beneath the KKNPP site (see Fig. 5a), two small aftershocks scenarios were simulated (see Table 5), whose hypocentres were located along the directrix of the Madonosaka syncline, at the two opposite sides with respect to its planar cross-section passing through the TEPCO facility (see Fig. 10a). Table 5. Summary of the aftershock parameters employed in this analysis. (ϕS; λ; δ) represents the strike, rake and dip angles estimated by F-NET Centroid Moment Tensor solution (Kubo et al.2002). Event  MJMA  M0[Nm]  (ϕS; λ; δ)[°]  τR[s]  AS1 (07/16/07-21:08)  4.4  52.1 · 1014  187; 70; 54  0.113  AS2 (07/16/07-17:42)  4.2  2.09 · 1014  309; 78; 37  0.045  Event  MJMA  M0[Nm]  (ϕS; λ; δ)[°]  τR[s]  AS1 (07/16/07-21:08)  4.4  52.1 · 1014  187; 70; 54  0.113  AS2 (07/16/07-17:42)  4.2  2.09 · 1014  309; 78; 37  0.045  View Large The aim is to characterize the 3-D effect of such a complex geology, by spanning different incidence angles of the wave front impinging the KKNPP. With respect to the LARGE model, the SMALL one, featured by a minimal shear-wave velocity of 700.0 m s−1 was designed to accurately propagate up to 5.0 Hz, with a third order polynomial degree. Except the eventual poor reproduction of the late surface wave arrivals, general features of observed records for AS1 and AS2 could be reproduced, as portrayed in Figs 11 (AS1, 1G1) and 13 (AS2, 5G1). Unfortunately, no recordings for AS1 at 5G1 are in the authors’ disposal. The improvement granted by the inclusion of the folded configuration however is evident (red traces in Figs 11b and d and 13b and d), compared to the poor fit to the records obtained with a subhorizontally layered geology not including the folding (green traces in Figs 11a and c, 12a and c and 13a and c). Figure 11. View largeDownload slide AS1: recorded (REC, blue) and synthetic velocigrams for subhorizontally layered not including the syncline-anticline structures (LAY, green) and folded (FLD, red) geology along the EW (top row) and NS (bottom row) directions, at 1G1. Figure 11. View largeDownload slide AS1: recorded (REC, blue) and synthetic velocigrams for subhorizontally layered not including the syncline-anticline structures (LAY, green) and folded (FLD, red) geology along the EW (top row) and NS (bottom row) directions, at 1G1. Figure 12. View largeDownload slide AS2: recorded (REC, blue) and synthetic velocigrams for subhorizontally layered not including the syncline-anticline structures (LAY, green) and folded (FLD, red) geology along the EW (top row) and NS (bottom row) directions, at 1G1. Figure 12. View largeDownload slide AS2: recorded (REC, blue) and synthetic velocigrams for subhorizontally layered not including the syncline-anticline structures (LAY, green) and folded (FLD, red) geology along the EW (top row) and NS (bottom row) directions, at 1G1. Figure 13. View largeDownload slide AS2: recorded (REC, blue) and synthetic velocigrams for subhorizontally layered not including the syncline-anticline structures (LAY, green) and folded (FLD, red) geology along the EW (top row) and NS (bottom row) directions, at 5G1. Figure 13. View largeDownload slide AS2: recorded (REC, blue) and synthetic velocigrams for subhorizontally layered not including the syncline-anticline structures (LAY, green) and folded (FLD, red) geology along the EW (top row) and NS (bottom row) directions, at 5G1. Figure 14. View largeDownload slide Scores of the Anderon’s Criteria along the EW/NS directions for AS1 (a and b) and AS2 (c and d). Figure 14. View largeDownload slide Scores of the Anderon’s Criteria along the EW/NS directions for AS1 (a and b) and AS2 (c and d). Figure 15. View largeDownload slide AS1: Comparison between recorded (REC, red for EW and orange for NS) response spectra Sa for SG4 and 1G1 sites compared to synthetic ones. (a and c) Sa spectra for layered geology not including the syncline-anticline structure; (b and d) Sa spectra for folded geology. Figure 15. View largeDownload slide AS1: Comparison between recorded (REC, red for EW and orange for NS) response spectra Sa for SG4 and 1G1 sites compared to synthetic ones. (a and c) Sa spectra for layered geology not including the syncline-anticline structure; (b and d) Sa spectra for folded geology. Figure 16. View largeDownload slide AS2: comparison between recorded (REC, red for EW and orange for NS) response spectra Sa for 5G1 and 1G1 sites compared to synthetic ones. (a and c) Sa spectra for layered geology not including the syncline-anticline structure; (b and d) Sa spectra for folded geology. Figure 16. View largeDownload slide AS2: comparison between recorded (REC, red for EW and orange for NS) response spectra Sa for 5G1 and 1G1 sites compared to synthetic ones. (a and c) Sa spectra for layered geology not including the syncline-anticline structure; (b and d) Sa spectra for folded geology. For both AS1 and AS2 and for both Unit 1 and Unit 5, the synthetic waveforms obtained with layered geology (LAY) at surface look out of phase with respect to the records, on one of the two directions EW/NS. This stress the great achievement obtained by including the Madonosaka syncline in the analysis. A more equitable judgment of the waveform fit is provided by the Anderson’s Criteria (Anderson 2004) in Figs 14(a) and (b) for AS1 at 1G1, and in Figs 14(c) and (d) for AS2 at 5G1. For both AS1 and AS2, the complex geological structure seems to influence principally the EW direction. However, the poor score obtained for both cases for the first four Anderson’s Criteria (related to the energy of the waveform) highlights one major drawback of the model: recordings are more energetic than synthetics, due to late arrivals that the model is not fit to reproduce. In terms of normalized Sa response spectra (with respect to their respective PGA), a suitable measure of the goodness of the model is the difference between the natural periods corresponding to the maximum Sa ordinate and referred as to TP. In Figs 15(b)–(d), the focusing of the wavefield due to the syncline causes a shift of TP towards the observed one, at both SG4 (Service Hall, G.L. −250 m) and 1G1 (Unit 1, G.L. 0 m), compared to the poor fit obtained for a subhorizontally layered geology not including the folding (see Figs 15a–c). For AS2, the pseudospectral acceleration response is less sensitive at Unit 5 (Figs 16a and b), whereas the fit is effectively improved at 1G1, by the presence of the folding in the model (Figs 16c and d). Finally, the KKNPP site response estimated by SEM3-D analyses for AS1 and AS2 is condensed in the following four figures: Figs 17 and 18 portray the pseudospectral acceleration spectra Sa estimated at KKNPP, for AS1 and AS2 respectively. Figure 17. View largeDownload slide AS1: Sa response spectra (in cm s−2) at different location around the KKNPP site. LAYERED (red) and FOLDED geological models were compared. Figure 17. View largeDownload slide AS1: Sa response spectra (in cm s−2) at different location around the KKNPP site. LAYERED (red) and FOLDED geological models were compared. Hayakawa et al. (2011) showed that the large amplitude of Unit 1 (located on the southwestward part of the KKNPP) comes from the folding structure where Unit 1 is located on the synclinal axis (Fig. 5a). In agreement with what proposed in the literature, Figs 17 and 18 indicate the folded geology as responsible of higher peaks for the response spectra obtained nearby Unit 1, for natural periods T < 0.5s. The site response is seemingly unaltered by the introduction of a complex folding structure beneath KKNPP (at Unit 5 and Service Hall). The minimum shear velocity introduced in the model corresponds to the engineering bedrock related to depth, justifying that the spatial incoherence of the ground motion recorded at surface and specifically the amplification southwestward are extremely influenced by the syncline–anticline structure, despite the effects due to shallow borehole geology (dispersion, attenuation/amplification) which have been disregarded in the analysis. This strengthens the argument that the site response at the Service Hall (which is located above Chuo-Yutai anticline) depended mainly on the nonlinear site-effects taking place at shallow depths (i.e. G.L.−250 m). The amplification trend occurs independently on the source position, being more accentuated along the EW direction. Watanabe et al. (2009) showed that the Upper Teradomari stratum does not alter the wave propagation of up-going waves, which tend to focalize at the Madonosaka syncline passing through Shiiya stratum. Site amplification becomes therefore significant at Unit 1 (KK1). On the other hand, Unit 5 (KK5, which is located above the Ushirodani anticline) is evidently more sensitive to wave-motion travelling from southwest to northeast, throughout the folding zone: the pseudospectral peaks differ from layered to folded geology only in Fig. 18, referring to AS2 (nucleated close by the third asperity on the fault rupture plane). Figure 18. View largeDownload slide Sa response spectra (in cm s−2) at different location around the KKNPP site. LAYERED (red) and FOLDED geological models were compared. Figure 18. View largeDownload slide Sa response spectra (in cm s−2) at different location around the KKNPP site. LAYERED (red) and FOLDED geological models were compared. 6 ON THE ACCURACY OF NUMERICAL SOLUTION The accuracy of the numerical method employed to solve the 3-D elastodynamic problem is notoriously a big deal to cope with. Due to inherent spatiotemporal nature of the wave-propagation physical phenomenon, the origins of the observed numerical dispersion are namely (1) the spatial and the (2) time discretization indeed. Usually, a Finite-Difference of Finite-Element spatial semi-discretization is assumed, along with a numerical time-integrator belonging to the Newmark’s family. In the SEM, high-order polynomials are employed, sampled at the GLL point on each element. Despite the intrinsic h − p property held by the SEM (leading to finer spatial discretizations and increased accuracy), it is intuitive that the two sources of numerical dispersions are not disconnected: Seriani & Oliveira (2008) stated that low-order time discretizations (such as the 2nd-order accurate Newmark’s schemes) deteriorate the high accuracy in space. According to Seriani & Oliveira (2008), P-wave dispersion and polarization errors are less sensitive to the value of Poisson’s ratio than the S-wave dispersion error. They confirmed that four grid points per wavelength are sufficient to have the dispersion error below 1 per cent on SE approximations of degree eight with GLL collocation points. Following the indications provided by the modal analysis performed by the mentioned authors, a parametric analysis on the computational model presented in this study for folded geology was performed. The minimum element size of the computational grid is ≈250 m at surface, with VS, min = 700.0 m s−1 and Poisson’s ratio ν ≈ 0.4. Three analysis were run therefore, with 5, 7 and 10 GLL points per element edge respectively. According to Seriani & Oliveira (2008) and for the Poisson’s ratio considered, the frequencies delimiting the accuracy of the analyses are 2.8, 5.0 and 7.0 Hz respectively. Therefore, the whole set of results of the three-analyses were low-pass filtered at 7 Hz, to check whether significant differences emerge. Fig. 19 shows the comparisons (in terms of accelerograms, Fourier’s spectrum along NS and EW directions respectively and geometric mean of the horizontal pseudospectral acceleration SaGMH) for 5, 7 and 10 GLL points respectively, at 1G1 (Unit 1, G.L. 0 m). The comparisons are rather explicative, in the following sense: the ground motion estimation does not change significantly for 7 or 10 GLL points, whereas it is rather underestimated (see the SaGMH ordinates in Fig. 19c) when the computational model is featured by 5 GLL points; the response at short periods is amplified by increasing the number of GLL points, granting a certain confidence to the analyses portrayed in the previous section (comparison with layered geology). Increasing the number of GLL points would lead to even more amplified response at both 1G1 and 5G1. the computational costs (measured in CPU-time, and plotted in the captioned subaxes in Fig. 19c) are growing exponentially with increasing number of GLL points, although no remarkable difference can be found between 7 and 10 GLL. This justify, for this particular case, the choice of 7 GLL as an optimum solution for the range of frequency of interest. Figure 19. View largeDownload slide Results of the parametric analysis on the SMALL model with folded geology at 1G1. (a) Accelerograms(NS direction), (b) Fourier’s spectra (NS direction) and (c) pseudospectral acceleration responses (geometric mean of horizontal components) for 5 (blue), 7 (red) and 10 (green) GLL respectively. Synthetics were filtered at 7 Hz. The computational costs of the SEM3D analyses are plotted in the captioned axes in (c). Figure 19. View largeDownload slide Results of the parametric analysis on the SMALL model with folded geology at 1G1. (a) Accelerograms(NS direction), (b) Fourier’s spectra (NS direction) and (c) pseudospectral acceleration responses (geometric mean of horizontal components) for 5 (blue), 7 (red) and 10 (green) GLL respectively. Synthetics were filtered at 7 Hz. The computational costs of the SEM3D analyses are plotted in the captioned axes in (c). 7 CONCLUSIONS In this paper, a source-to-site numerical model of the NCOEQ-2007 earthquake is constructed, to reproduce the seismic site response observed at KKNPP. Specifically, the aim of this study is to criticize the common simplified assumption of a subhorizontally layered Earth’s crust, inferred from seismic and geodetic and geological surveys, stressing the importance of the shallow geological conformation in altering the incident wave motion. The case of KKNPP rises as an exemplary benchmark, due to extensively observed near-field effects at the site, such as the ground motion incoherence within a relatively small distance. As a matter of fact, the incident wave motion resulted amplified in the south-west area, nearby the first group of nuclear reactors and structures (Unit 1–4), compared to the second group of facilities (Unit 5–7), located northeastward. Two aftershocks belonging to the NCOEQ-2007 sequence were simulated to prove the assumption focusing on the effect of the geological configuration. The dense observation network deployed in the Niigata region and at the KKNPP steered the calibration of two numerical seismic scenarios: one to portrays the regional incident wavefield (in a frequency band 0.1–0.5 Hz) and a second one restricted to the ground shaking pertaining the KKNPP surroundings, but extending the modelling accuracy up to 5.0 Hz (and possibly 7.0 Hz). The paper explains the methodological approach to construct an holistic source-to-site computational model: LARGE model is first calibrated upon the results of semi-analytical solutions, considering a 1-D regional geology so to constrain the low-frequency band of the incident wavefield; SMALL model is later constructed by trimming a portion of LARGE model and plugging into the 1-D regional geology the intricate 3-D set of syncline-anticlines, placed right underneath KKNPP. LARGE model is used to assess the effect of the source parameters and of the topography, comparing the results with classical GMPEs. SMALL model includes the Japan Sea, the bathymetry and the coastline, other than the folding geology. The outcome of the numerical exercise provides a fairly good agreement between the recordings and the synthetics at KKNPP, with high accuracy at 5.0 Hz. The effect of the folding clearly improves this fit, compared to the too much simplified subhorizontally layered profiles. Moreover, the effort of including a complicated geology in the seismological model is worthwhile, as already proved by Tsuda et al. (2011), but herein presented in a more extended way, including rigorous comparison between time-histories, goodness of fit estimation at several locations and detailed seismic site response at KKNPP. One major contribution of this paper is that the two-step construction of a hybrid velocity structure bridges the gap between previous studies targeting the regional wavefield (but limited in frequency, e.g. Aochi et al.2013) and more site-specific analyses, describing the wavefield at KKNPP (Watanabe et al.2009; Tsuda et al.2011). The numerical analyses highlight the impact of the syncline laying below Unit 1 in causing an amplification (both in terms of peak ground motions and response spectra) with respect to the layered geology. The syncline seemingly focalizes the upgoing wavefield, drifting the radiated energy towards Unit 1. This occurrence represents an exception to the current design standards, which traditionally leans towards the assumption of a 1-D layered geology in case of lack of further information concerning the subsurface geology. The case study stresses the importance of the propagation path in drifting and altering (eventually amplifying) the wavefield radiated even from relatively weak aftershocks. Those findings confirm previous studies and extend their validity at higher frequency (5.0 Hz). Moreover, in this paper a major effort of quantifying the numerical accuracy is performed and presented aside the main analysis, completing the results and depicting the uncertainty related. However, the major shortcoming of the analyses presented herein resides in the lack of more refined models of the scattered wavefield in a high frequency range (0–10 Hz). This is a crucial feature to be taken into consideration in the transition towards a broad-band deterministic modelling of the earthquake phenomenon. For instance, the inclusion of heterogeneous soil deposits and nonlinear hysteretic damping would surely enhance the realism of the final outcome, although it may mislead the understanding of the physical mechanisms taking place, due to the increased number of parameters required. Based on these considerations, it can be argued that the current seismic design guidelines for critical structures should require further improvements of site-characterization not only in terms of active faults and local geotechnical conditions, but specifically on the conformation of the Earth’s crust (such as the SCEC research plan for the period 2012-2016, focused on specific areas in the Southern California and enhanced the Community Velocity Model (CVM), describing seismic P- and S-wave velocities and densities throughout the southern California region). It is in the authors’ belief that the modern procedures for earthquake ground motion prediction for seismic hazard analysis, especially for the low probability events inherent to the seismic design of critical structures, should account for a reliable incident wavefield issued from forward deterministic analyses. The ever increasing computational power leads to put major efforts to extend the frequency band and to enable PBS to be used with confidence in engineering applications (e.g. seismic hazard analysis and structural or geotechnical dynamic analysis). Acknowledgements This work, within the SINAPS@ project, benefited from French state funding managed by the National Research Agency under program RNSR Future Investments bearing reference No. ANR-11-RSNR-0022-04. The research reported in this paper has been supported in part by the SEISM Paris Saclay Research Institute. The authors are very grateful to the Tokyo Electric Power Company (TEPCO) for providing such high-quality earthquake recordings of the Niigata-Ken Chūetsu-Oki, reported in TEPCO (2007). Time histories and velocity profiles used in this study were collected from the KiK-net website: http://www.kik.bosai.go.jp/kik/ (last accessed November 2011). The authors are very grateful to the National Research Institute for Earth Science and Disaster Prevention (NIED) for providing such high-quality earthquake recordings. The authors are very thankful to the Data Management Center of the National Research Institute for Earth Science and Disaster Resilience (NIED DMC), which provides continuous seismic waveform and other seismological data obtained by the NIED High Sensitivity Seismograph Network (NIED Hi-net) and the other seismograph networks in Japan. NIED, Japan Meteorological Agency (JMA), and universities have started exchanging high-sensitivity seismograph data obtained by each seismic network with a specific policy for ‘Exchange and Open Access of Waveform Data obtained by High Sensitivity Seismograph Observation’ conducted by the Headquarters for Earthquake Research Promotion of the Japanese government. The website (http://www.hinet.bosai.go.jp/?LANG = en) provides not only the NIED Hi-net waveform data but also the other high-sensitivity seismograph data. The Digital Elevation Model belongs to the Shuttle Radar Topography Mission (SRTM) database, result from a collaborative effort by the National Aeronautics and Space Administration (NASA) and the National Geospatial-Intelligence Agency (NGA), as well as the participation of the German and Italian space agencies, to generate a near-global digital elevation model (DEM) of the Earth using radar interferometry. Source: https://dds.cr.usgs.gov/srtm/version2_1. REFERENCES Aki K., 1967. Scaling law of seismic spectrum, J. geophys. Res. , 72( 4), 1217– 1231. https://doi.org/10.1029/JZ072i004p01217 Google Scholar CrossRef Search ADS   Anderson J.G., 2004. Quantitative measure of the goodness-of-fit of synthetic seismograms, in13th World Conference on Earthquake Engineering, Vancouver, BC, Canada, Paper No. 243, pp. 775– 784. Aochi H., Yoshimi M., 2016. Seismological asperities from the point of view of dynamic rupture modeling: the 2007 Mw6.6 Chuetsu-Oki, Japan, earthquake, J. Seismol. , 20( 4), 1089– 1105. https://doi.org/10.1007/s10950-016-9569-5 Google Scholar CrossRef Search ADS PubMed  Aochi H., Ducellier A., Dupros F., Delatre M., Ulrich T., de Martin F., Yoshimi M., 2013. Finite difference simulations of seismic wave propagation for the 2007 Mw6.6 Niigata-ken Chuetsu-Oki earthquake: validity of models and reliable input ground motion in the near field, Pure appl. Geophys. , 170( 1–2), 43– 64. https://doi.org/10.1007/s00024-011-0429-5 Google Scholar CrossRef Search ADS   Baker J.W., 2007. Quantitative classification of near-fault ground motions using wavelet analysis, Bull. seism. Soc. Am. , 97( 5), 1486– 1501. https://doi.org/10.1785/0120060255 Google Scholar CrossRef Search ADS   Bao H., Bielak J., Ghattas O., Kallivokas L., O’Hallaron D., Shewchuk J., Xu J., 1998. Large-scale simulation of elastic wave propagation in heterogeneous media on parallel computers, Comput. Methods Appl. Mech. Eng. , 152( 1), 85– 102. https://doi.org/10.1016/S0045-7825(97)00183-7 Google Scholar CrossRef Search ADS   Bard P.-Y., Bouchon M., 1980. The seismic response of sediment-filled valleys. Part 1. The case of incident SH waves, Bull. seism. Soc. Am. , 70( 4), 1263– 1286. Bizzarri A., 2014. On the point-source approximation of earthquake dynamics, Ann. Geophys. , 57( 3). doi: http://dx.doi.org/10.4401/ag-6479. Bouchon M., 1981. A simple method to calculate Green’s functions for elastic layered media, Bull. seism. Soc. Am. , 71( 4), 959– 971. Brune J.N., 1970. Tectonic stress and the spectra of seismic shear waves from earthquakes, J. geophys. Res. , 75( 26), 4997– 5009. https://doi.org/10.1029/JB075i026p04997 Google Scholar CrossRef Search ADS   Chaljub E. et al.  , 2015. 3-D numerical simulations of earthquake ground motion in sedimentary basins: testing accuracy through stringent models, Geophys. J. Int. , 201( 1), 90– 111. https://doi.org/10.1093/gji/ggu472 Google Scholar CrossRef Search ADS   Cirella A., Piatanesi A., Tinti E., Cocco M., 2008. Rupture process of the 2007 Niigataken Chuetsu-Oki earthquake by non-linear joint inversion of strong motion and GPS data, Geophys. Res. Lett. , 35( 16), 1– 5. https://doi.org/10.1029/2008GL034756 Google Scholar CrossRef Search ADS   Cotton F., Pousse G., Bonilla F., Scherbaum F., 2008. On the discrepancy of recent European ground-motion observations and predictions from empirical models: Analysis of KiK-net accelerometric data and point-sources stochastic simulations, Bull. seism. Soc. Am. , 98( 5), 2244– 2261. https://doi.org/10.1785/0120060084 Google Scholar CrossRef Search ADS   Cox K.E., Ashford S.A., 2002. Characterization of large velocity pulses for laboratory testing, Tech. Rep. April , Pacific Earthquake Engineering Research Center. Day S.M., Graves R., Bielak J., Dreger D., Larsen S., Olsen K.B., Pitarka A., Ramirez-Guzman L., 2008. Model for basin effects on long-period response spectra in southern California, Earthq. Spectra , 24( 1), 257– 277. https://doi.org/10.1193/1.2857545 Google Scholar CrossRef Search ADS   De Martin F., 2011. Verification of a spectral-element method code for the Southern California earthquake center LOH.3 viscoelastic case, Bull. seism. Soc. Am. , 101( 6), 2855– 2865. https://doi.org/10.1785/0120100305 Google Scholar CrossRef Search ADS   De Martin F., Aochi H., Foerster E., 2007. Testing the Double-Couple Point Source Model in GEFDYN, Tech. rep. , BRGM. Dreger D., Hurtado G., Chopra A.K., Larsen S., 2007. Near-Fault Seismic Ground Motions, Tech. Rep. 59A0435 , California Department of Transportation. Ducellier A., Aochi H., 2010. Numerical simulation of the Mw6.6 Niigata, Japan, earthqhake: reliable input ground motion for engineering purpose, in 14th European Conference on Earthquake Engineering , Ohrid, Macedonia. Dupros F., De Martin F., Foerster E., Komatitsch D., Roman J., 2010. High-performance finite-element simulations of seismic wave propagation in three-dimensional nonlinear inelastic geological media, Parallel Comput. , 36( 5-6), 308– 325. https://doi.org/10.1016/j.parco.2009.12.011 Google Scholar CrossRef Search ADS   Duputel Z., Rivera L., Kanamori H., Hayes G., 2012. W phase source inversion for moderate to large earthquakes (1990-2010), Geophys. J. Int. , 189( 2), 1125. https://doi.org/10.1111/j.1365-246X.2012.05419.x Google Scholar CrossRef Search ADS   Duputel Z., Tsai V.C., Rivera L., Kanamori H., 2013. Using centroid time-delays to characterize source durations and identify earthquakes with unique characteristics, Earth planet. Sci. Lett. , 374, 92– 100. https://doi.org/10.1016/j.epsl.2013.05.024 Google Scholar CrossRef Search ADS   Faccioli E., Vanini M., 2003. Complex seismic site effects in sediment-filled valleys and implications on design spectra, Prog. Struct. Eng. Mater. , 5( 4), 223– 238. https://doi.org/10.1002/pse.156 Google Scholar CrossRef Search ADS   Faccioli E., Maggio F., Paolucci R., Quarteroni A., 1997. 2D and 3D elastic wave propagation by a pseudo-spectral domain decomposition method, J. Seismol. , 1( 3), 237– 251. https://doi.org/10.1023/A:1009758820546 Google Scholar CrossRef Search ADS   Fujiwara H. et al.  , 2009. A study on subsurface structure model for deep sedimentary layers of Japan for strong-motion evaluation, Technical Note of the National Research Institute for Earth Science and Disaster Prevention , No. 337, 260. Fukushima S., Hayashi T., Yashiro H., 2007. Seismic hazard analysis based on the joint probability density function of PGA and PGV, in Transactions, SMiRT 19 , Vol. M03. Gallovič F., Käser M., Burjánek J., Papaioannou C., 2010. Three-dimensional modeling of near-fault ground motions with nonplanar rupture models and topography: case of the 2004 Parkfield earthquake, J. geophys. Res. , 115( B3), 1– 17. Google Scholar CrossRef Search ADS   Gélis C., Bonilla L.F., 2012. 2-D P-SV numerical study of soil-source interaction in a non-linear basin, Geophys. J. Int. , 191, 1374– 1390. Göddeke D., Komatitsch D., Möller M., 2014. Finite and spectral element methods on unstructured grids for flow and wave propagation methods, chap.9, in Numerical Computations with GPUs , pp. 183– 206, ed. Kindratenko V., Springer. Graves R.W., 1996. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences, Bull. seism. Soc. Am. , 86( 4), 1091– 1106. Graves R.W., Wald D.J., 2004. Observed and simulated ground motions in the San Bernardino Basin Region for the Hector Mine, California, earthquake, Bull. seism. Soc. Am. , 94( 1), 131– 146. https://doi.org/10.1785/0120030025 Google Scholar CrossRef Search ADS   Gürpinar A., Serva L., Livio F., Rizzo P.C., 2017. Earthquake-induced crustal deformation and consequences for fault displacement hazard analysis of nuclear power plants, Nucl. Eng. Des. , 311, 69– 85. https://doi.org/10.1016/j.nucengdes.2016.11.007 Google Scholar CrossRef Search ADS   Hartzell S.H., 1978. Earthquake aftershocks as Green’s functions, Geophys. Res. Lett. , 5( 1), 1– 4. https://doi.org/10.1029/GL005i001p00001 Google Scholar CrossRef Search ADS   Hayakawa T., Tsuda K., Uetake T., Hikima K., Tokumitsu R., Nagumo H., 2011. Modeling 3D velocity structure in the fault region of the 2007 Niigataken Chuetu-oki Earthquake-Incorporating the 3D fold geological structure beneath the Kashiwazaki-Kariwa nuclear power plant, in Japan Geoscience Union Meeting , no. SSS023, p. 14. Hijikata K., Nishimura I.and Mizutani H., Tokumitsu R., Mashimo M., Tanaka S., 2010. Ground motion characteristics of 2007 Niigata-ken Chuetsu-Oki earthquake, J. Struct. Construct. Eng. , 75( 653), 1279– 1288. https://doi.org/10.3130/aijs.75.1279 Google Scholar CrossRef Search ADS   Hisada Y., 1994. An efficient method for computing Green’s Functions for layered half-space with sources and receivers at close depths, Bull. seism. Soc. Am. , 84( 5), 1456– 1472. Hisada Y., 1995. An efficient method for computing Green’s functions for a layered half-space with sources and receivers at close depths (part 2), Bull. seism. Soc. Am. , 85( 4), 1080– 1093. Hisada Y., 2008. Broadband strong motion simulation in layered half-space using stochastic Green’s function technique, J. Seismol. , 12( 2), 265– 279. https://doi.org/10.1007/s10950-008-9090-6 Google Scholar CrossRef Search ADS   Irikura K., 1983. Semi-Empirical Estimation of Strong Ground Motions During, Bull. Disaster Prevention Res. Inst. , 33( 2), 63– 104. JNES, 2008. Construction of the subsurface structure around the source area of 2007 Niigata Ken Chuetsu-Oki earthquake, Tech. rep. , Japanes Nuclear Energy Safety Organization. Kato A., Sakai S., Kurashimo E., Igarashi T., Iidaka T., Hirata N., Iwasaki T., Kanazawa T., 2008. Imaging heterogeneous velocity structures and complex aftershock distributions in the source region of the 2007 Niigataken Chuetsu-Oki Earthquake by a dense seismic observation, Earth Planets Space , 60, 1111– 1116. https://doi.org/10.1186/BF03353145 Google Scholar CrossRef Search ADS   Kato A. et al.  , 2009. Reactivation of ancient rift systems triggers devastating intraplate earthquakes, Geophys. Res. Lett. , 36( L05301), 1– 5. Kawabe H., Kamae K., 2010. Source modelling and 3D ground motion simulation of the 2007 Niigataken Chuetsu-oki earthquake (Mj6.8), in 13th Japan Earthquake Engineering Symposium , Tsukuba Japan. Kayen R. et al.  , 2007. Investigation of the M6.6 Niigata-Chuetsu Oki, Japan, Earthquake of July 16, 2007, Tech. Rep. 1365 , U.S. Department of the Interior and U.S. Geological Survey (USGS). Kobayashi I., Tateishi M., Yoshimura T., Ueda T., Kato H., 1995. Geology of the Kashiwazakii district. Quadrangle Series: Scale 1:50000, Vol. 7, Issue 37, P. 101, Geological Survey of Japan. Komatitsch D., Vilotte J.-P., 1998. The spectral element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures, Bull. seism. Soc. Am. , 88( 2), 368– 392. Komatitsch D., Vilotte J.P., Vai R., Castillo-Covarrubias J.M., Sánchez-Sesma F.J., 1999. The spectral element method for elastic wave equations: application to 2D and 3D seismic problems, in SEG Expanded Abstracts 17 , vol. 45, pp. 1139– 1164. Kubo A., Fukuyama E., Kawai H., Nonomura K., 2002. NIED seismic moment tensor catalogue for regional earthquakes around Japan: quality test and application, Tectonophysics , 356( 1), 23– 48. https://doi.org/10.1016/S0040-1951(02)00375-X Google Scholar CrossRef Search ADS   Mazzieri I., Smerzini C., Antonietti P.F., Rapetti F., Stupazzini M., Paolucci R., Quarteroni A., 2011. Non-conforming spectral approximations for the elastic wave equation in heterogeneous media, in COMPDYN 2011 ECCOMAS Thematic Conference on Computational Methods in Structural Dynamics and Earthquake Engineering , pp. 1– 17, eds Papadrakakis M., Fragiadakis M., Plevris V., Corfu, Greece. Mazzieri I., Stupazzini M., Guidotti R., Smerzini C., 2013. SPEED: SPectral Elements in Elastodynamics with Discontinuous Galerkin: a non-conforming approach for 3D multi-scale problems, Int. J. Numer. Methods Eng. , 95( 12), 991– 1010. https://doi.org/10.1002/nme.4532 Google Scholar CrossRef Search ADS   Miyake H., Koketsu K., Hikima K., Shinohara M., Kanazawa T., 2010. Source Fault of the 2007 Chuetsu-Oki, Japan, Earthquake, Bull. seism. Soc. Am. , 100( 1), 384– 391. https://doi.org/10.1785/0120090126 Google Scholar CrossRef Search ADS   Moczo P., Kristek J., Galis M., Pazak P., Balazovjech M., 2007. The finite-difference and finite-element modeling of seismic wave propagation and earthquake motion, Acta Phys. Slovaca , 57( 2), 177– 406. https://doi.org/10.2478/v10155-010-0084-x Google Scholar CrossRef Search ADS   Obara K., Kasahara K., Hori S., Okada Y., 2005. A densely distributed high-sensitivity seismograph network in Japan: Hi-net by National Research Institute for Earth Science and Disaster Prevention, Rev. Sci. Instrum. , 76( 2), 1– 12. https://doi.org/10.1063/1.1854197 Google Scholar CrossRef Search ADS   Paolucci R., Mazzieri I., Smerzini C., Stupazzini M., 2014. Physics -Based Earthquake Ground Shaking Scenarios in Large Urban Areas, in Perspectives on European Earthquake Engineering and Seismology , Vol. 34 of Geotechnical, Geological and Earthquake Engineering, pp. 331– 359, ed. Ansal A., Springer. Google Scholar CrossRef Search ADS   Paolucci R., Mazzieri I., Smerzini C., 2015. Anatomy of strong ground motion: near-source records and 3D physics-based numerical simulations of the Mw 6.0 May 29 2012 Po Plain earthquake, Italy, Geophys. J. Int. , 203, 2001– 2020. https://doi.org/10.1093/gji/ggv405 Google Scholar CrossRef Search ADS   Pavlenko O.V., Irikura K., 2012. Nonlinear soil behavior at the Kashiwazaki-Kariwa Nuclear Power Plant during the Niigata Chuetsu-Oki earthquake (July, 16, 2007), Pure appl. Geophys. , 169( 10), 1777– 1800. https://doi.org/10.1007/s00024-011-0447-3 Google Scholar CrossRef Search ADS   Pitarka A., Irikura K., Iwata T., Sekiguchi H., 1998. Three-dimensional simulation of the near-fault ground motion for the 1995 Hyogo-Ken Nanbu (Kobe), Japan, earthquake, Bull. seism. Soc. Am. , 88( 2), 428– 440. Quinay P. E.B., Ichimura T., Hori M., Nishida A., Yoshimura S., 2013. Seismic Structural Response Estimates of a Fault-Structure System Model with Fine Resolution Using Multiscale Analysis with Parallel Simulation of Seismic-Wave Propagation, Bull. seism. Soc. Am. , 103( 3), 2094– 2110. https://doi.org/10.1785/0120120216 Google Scholar CrossRef Search ADS   Sekiguchi H. et al.  , 2009. 3D subsurface structure model of the Niigata sedimentary basin, Ann. Rep. Active Fault Paleoearthq. Res. , 9, 175– 259. Seriani G., Oliveira S.P., 2008. Dispersion analysis of spectral element methods for elastic wave propagation, Wave Motion , 45( 6), 729– 744. https://doi.org/10.1016/j.wavemoti.2007.11.007 Google Scholar CrossRef Search ADS   Shiba Y., 2008. Source process and broadband strong motions during the Niigata-ken Chuetsu-Oki earthquake in 2007, Denryoku Chuo Kenkyusho Hokoku , ( N08007), 1– 4, Japanese with English Abstract. Shinohara M. et al.  , 2008. Precise aftershock distribution of the 2007 Chuetsu-Oki Earthquake obtained by using an ocean bottom seismometer network, Earth Planets Space , 60( 11), 1121– 1126. https://doi.org/10.1186/BF03353147 Google Scholar CrossRef Search ADS   Smerzini C., Villani M., 2012. Broadband Numerical Simulations in Complex Near-Field Geological Configurations: The Case of the 2009 Mw 6.3 L’Aquila Earthquake, Bull. seism. Soc. Am. , 102( 6), 2436– 2451. https://doi.org/10.1785/0120120002 Google Scholar CrossRef Search ADS   Spudich P., Archuleta R., 1987. Techniques for earthquake ground-motion calculation with applications to source parameterization to finite faults, in Seismic Strong Motion Synthesis , chap. pp. 205– 265, ed. Bolt B.A., Academic Press. Google Scholar CrossRef Search ADS   Taborda R., Bielak J., 2011. Large-scale earthquake simulation: computational seismology and complex engineering systems, Comput. Sci. Eng. , 13( 4), 14– 27. https://doi.org/10.1109/MCSE.2011.19 Google Scholar CrossRef Search ADS   Taborda R., Bielak J., 2013. Ground-motion simulation and validation of the 2008 Chino Hills, California, earthquake, Bull. seism. Soc. Am. , 103( 1), 131– 156. https://doi.org/10.1785/0120110325 Google Scholar CrossRef Search ADS   TEPCO, 2007. The data analysis recorded at the Kashiwazaki Kariwa Nuclear Power Plant during the 2007 Niigata-ken Chuetsu-oki earthquake, Tech. rep. , The Tokyo Electric Power Company, Inc, in Japanese. Tokumitsu R., Kikuchi M., Nishimura I., Shiba Y., Tanaka S., 2009. Analysis of the strong motion records obtained from the 2007 Niigataken Chuetsu-Oki earthquake and determination of the design basis ground motions at the Kashiwazaki Kariwa Nuclear Power Plant. Part 1. Outline of the strong motion records and estimation of factors in large amplification, Tech. rep. Tsuda K., Hayakawa T., Uetake T., Hikima K., Tokimitsu R., Nagumo H., Shiba Y., 2011. Modeling 3D Velocity Structure in the Fault Region of the 2007 Niigataken Chuetu-Oki Earthquake with Folding Structure, in 4th IASPEI/IAEE International Symposium on Effects of Surface Geology on Seismic Motion , Santa Barbara, CA, pp. 1– 11. Uetake T., Nishimura I., Mizutani H., 2008. Characteristics of strong motion records in Kashiwazaki-Kariwa Nuclear Power Station during the Niigataken Chuetsu-Oki earthquake in 2007, in 14th WCEE , Beijing, China. Watanabe T., Moroi T., Nagano M., Tokumitsu R., Kikuchi M., Nishimura I., 2009. Analysis of the strong motion records obtained from the 2007 Niigataken Chuetsu-Oki earthquake and determination of the design basis ground motions at the Kashiwazaki Kariwa Nuclear Power Plant. Part 2. Difference of site amplification based on the 2D FEM analysis of the folded structure, Tech. rep. Watanabe T., Moroi T., Tokumitsu R., Nishimura I., Hijikata K., 2011. Examination of relation between locations of asperities and site amplification characteristics of ground motions by analysis considering the folded structure-Estimation based on the strong motion records obtained from the 2007 Niigataken Chuetsu-oki earthquake in the Kashiwazaki-Kariwa Nuclear Power Station, J. Struct. Construct. Eng. , 76( 659), 71– 78. https://doi.org/10.3130/aijs.76.71 Google Scholar CrossRef Search ADS   Yee E., Stewart J., Tokimatsu K., 2011. Nonlinear Site Response and Seismic Compression at Vertical Array Strongly Shaken by the Niigata-ken Chuetsu-Oki Earthquake, Tech. rep. , Pacific Earthquake Engineering Research Center-College of Engineering, University of California, Berkeley. APPENDIX: To assess the goodness of fit, the Anderson’s criteria are employed (Anderson 2004, listed in Table A1). Table A1. Goodness of fit parameters. Number  Symbol  Similarity of:  C1  SDa  Arias duration  C2  SDe  Energy duration  C3  SIa  Arias Intensity  C4  SIv  Energy Integral  C5  Spga  Peak Acceleration  C6  Spgv  Peak Velocity  C7  Spgd  Peak Displacement  C8  Ssa  Response Spectra  C9  Sfs  Fourier Spectra  C10  C*  Cross Correlation  Number  Symbol  Similarity of:  C1  SDa  Arias duration  C2  SDe  Energy duration  C3  SIa  Arias Intensity  C4  SIv  Energy Integral  C5  Spga  Peak Acceleration  C6  Spgv  Peak Velocity  C7  Spgd  Peak Displacement  C8  Ssa  Response Spectra  C9  Sfs  Fourier Spectra  C10  C*  Cross Correlation  View Large © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations