David Castro-Cruz

Filippo Gatti

Fernando Lopez-Caballero

Assessing the impact of regional geology on the ground motion model variability at the Kashiwazaki-Kariwa Nuclear Power Plant (Japan) via physics-based numerical simulation

Keywords:

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

With the increasing broader access to high-quality seismic databases, empirical relationships to predict ground motion intensity flourished, cast into Ground Motion Prediction Equations (GMPEs, whose omni-comprehensive compendium can be found at http://www.gmpe.org.uk/, see [START_REF] Douglas | Ground motion prediction equations 1964-2021[END_REF], of the general form:

log = { } ; { } + Δ (1)
where is the generic intensity measure (IM) for earthquake at station . The pseudo-spectral acceleration (approximating the maximum absolute acceleration of the single degree of freedom system with 5% damping), the Peak Ground Acceleration (PGA) and the Peak Ground Velocity (PGV) are among the most commonly used IMs. The functional form embodies the geometric mean IM prediction (log value). { } embraces the source characteristics (usually earthquake magnitude, type of earthquake, focal mechanism), the notion of source-to-site distance (hypocentral, epicentral and Joyner-Boore among others), the site parameters (usually the site class based on the average shear-wave velocity of the first 30 m below the ground surface, noted as V ,30 ). { } is the set of coefficients obtained via non-linear regression on large seismic databases. The residual Δ represents the shift between the geometric mean prediction and observations. Δ is usually modeled as a centered normal random variable with distribution  (0, ) [START_REF] Atik | The Variability of Ground-Motion Prediction Models and Its Components[END_REF][START_REF] Rodriguez-Marek | A Model for Single-Station Standard Deviation Using Data from Various Tectonic Regions[END_REF]. Moreover, Δ is modeled as the sum of the between-event variability (the earthquake-wise station-average deviation of the observed ground motion from the geometric mean ) and of the within-event variability

(the stationwise misfit between the observed earthquake ground motion and the sum + , [START_REF] Rodriguez-Marek | A Model for Single-Station Standard Deviation Using Data from Various Tectonic Regions[END_REF].

Δ = + (2)
and are both modeled as centered normal random variable with distributions  (0, ) and  (0, ) respectively.

GMPEs usually benefit from the ever-increasing access to seismic data and meta-data. More complex functional can be calibrated by adding new parameters to the recipe in order to take into account site-specific and source-specific features. However, current seismic hazard analysis struggles in reducing the uncertainty Δ and the associated standard deviation (especially for earthquakes with high return periods, [START_REF] Strasser | Sigma: Issues, Insights, and Challenges[END_REF]. Some authors achieved this goal by removing the so-called ergodicity assumption (see, for instance, [START_REF] Atik | The Variability of Ground-Motion Prediction Models and Its Components[END_REF], i.e. including site-and scenario-dependent parameters. However, this approach suffers from (1) a lack of observed records in some regions to sufficiently constrain the regression process and (2) poor generalization, i.e., the impossibility of extrapolating reliable predictions in general contexts [START_REF] Abrahamson | Application of single-station sigma ground motion prediction equations in practice[END_REF][START_REF] Chen | A new method for estimation of the attenuation relationship with variance components[END_REF]. Those drawbacks hinder the establishment of high-fidelity site-/scenario-specific GMPEs. A very appealing alternative is provided by high-fidelity earthquake simulations, duly validated against past recorded earthquakes [START_REF] Bradley | On-going challenges in physics-based ground motion prediction and insights from the 2010-2011 Canterbury and 2016 Kaikoura, New Zealand earthquakes[END_REF]. This approach sacrifices the GMPE handy regression framework to explicitly tackle the modeling uncertainties at stake, namely: (1)

the definition of a reasonable and realistic rupture scenario on active faults is required, (2) the 3-D crustal structure, (3) the surficial site conditions (related to the so-called site effects). For those reasons, waveform inversion and array backprojection are employed to detect and characterize fault asperities that generated the strong ground motion [START_REF] Aoi | Source process of the 2007 Niigata-ken Chuetsu-oki earthquake derived from nearfault strong motion data[END_REF][START_REF] Honda | Array back-projection imaging of the 2007 Niigataken Chuetsu-Oki earthquake striking the world's largest nuclear power plant[END_REF]. GMPEs and numerical simulations are being used together for Probabilistic Seismic Hazard Assessment (PSHA). For instance, [START_REF] Milner | Toward PhysicsâĂŘBased Nonergodic PSHA: A Prototype Fully Deterministic Seismic Hazard Model for Southern California[END_REF] managed by updating non-ergodic PSHA for Southern California via direct numerical simulation.

The identification of the geology in the region of interest (bedrock and softer sediments) is either achieved by interpretation and analysis of borehole logs in the surroundings [START_REF] Fantoni | Tectono-sedimentary setting of the Po Plain and Adriatic foreland[END_REF], or by seismic tomography [START_REF] Lee | Full 3D tomography for crustal structure in Southern California based on the scattering integral and the adjoint wavefield methods[END_REF][START_REF] Jian | Prediction of oil-bearing single sandbody by 3d geological modeling combined with seismic inversion[END_REF]. Finally, laboratory tests are necessary when including non-linear soil rheology (see the PRENOLIN international benchmark on 1-D numerical wave-propagation in [START_REF] Régnier | International Benchmark on Numerical Simulations for 1D, Nonlinear Site Response (PRENOLIN): Verification Phase Based on Canonical Cases[END_REF]. Despite the computational performance achieved by the numerical tools, the mentioned uncertainties prevent the complete characterization of broad-band ground shaking scenarios, especially at short periods (i.e. < 0.5 s, [START_REF] Paolucci | Broad-band ground motions from 3D physics-based numerical simulations using Artificial Neural Networks[END_REF]. In recent years, physics-based numerical simulation improved by merging its fidelity at long periods (usually for natural periods > 0.2 s) with data-driven ground motion models, conditioned by the outcome of the numerical analysis [START_REF] Gatti | Towards blending Physics-Based numerical simulations and seismic databases using Generative Adversarial Network[END_REF][START_REF] Castro-Cruz | High-fidelity broad-band prediction of regional seismic response: a hybrid coupling of Physics-Based synthetic simulation and Empirical Green's functions[END_REF][START_REF] Jayalakshmi | Hybrid broadband ground motion simulations in the indo-gangetic basin for great himalayan earthquake scenarios[END_REF]. Several comparisons [START_REF] Olsen | 3D Viscoelastic Wave Propagation in the Upper Borrego Valley, California, Constrained by Borehole and Surface Data[END_REF][START_REF] Smerzini | Comparison of 3D, 2D and 1D numerical approaches to predict long period earthquake ground motion in the Gubbio plain, Central Italy[END_REF] showed large differences between 1-D 2-D and 3-D simulations, establishing that 3-D simulations can better capture the source footprint and the impact of the geological structures. However, the power of 3-D models is sometimes limited by the only available 1-D velocity profiles around specific areas.

In this work, high-fidelity 3-D numerical simulations are exploited to clarify their sensitivity to the geological model considered in a relatively low-frequency range (0-3 Hz). The uncertainty quantification is inferred from the reproduced variability of the ground motions for three different yet plausible geological models. The analysis follows the context of GMPE predictions, given the large number of synthetic time-histories extracted from the numerical simulation at several geographical locations in the region of interest. This work focuses on a real case study, the Niigata prefecture (central-west Japan), which hosts the largest nuclear power plant in the world, located in Kashiwazaki-Kariwa (KKNPP, see Figure 1).

This region has been struck by several significant earthquakes, with the remarkable 6.6 Niigata Chuetsu-Oki earthquake of 2007 (NCOEQ2007), which did cause the KKNPP shut-down. Gatti et al. (2018a) showed the importance of considering the complex 3-D shallow underground structure to approach the recorded seismic response, compared to traditional layered geology. They confirmed an improved fit of the recorded time-histories at KKNPP, when the complex folded geology below KKNPP was included in their numerical model. Gatti et al. (2018c) employed the technique named ANN2BB [START_REF] Paolucci | Broad-band ground motions from 3D physics-based numerical simulations using Artificial Neural Networks[END_REF] to enrich this synthetic wave-motion at high-frequency and successfully perform a Soil-Structure Interaction (SSI) study of the Unit 7 reactor building at KKNPP. These studies highlighted the difficulty of constructing and calibrating a regional digital twin, including the uncertainty of the available geological information (usually 2-D geological maps) and the uncertainty related to the numerical discretization.

Data

3-D numerical models require detailed information on the site geological configuration and source characteristics, in order to accurately reproduce the complex 3-D wave field. Ideally, high-quality information is required at three different scales: the regional scale (strictly connected to the long-period part of the rendered ground motion), the site-scale (at which most of the non-linear multi-path effects take place), and the borehole scale (necessary when investigating complex SSI effects). Nevertheless, these ideal conditions are rarely met, given the tremendous effort needed to collect seismological observations, to perform geophysical and geological campaigns, to invert and to exploit those data. More realistically, 1-D geological profiles are provided for the seismological region (characteristic dimension ≈ 1 km). 2-D or 3-D geological structures are also available for some sites of interest, such as at the KKNPP [START_REF] Watanabe | 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[END_REF].

Available 1-D velocity structures

Figure 2 compares the available 1-D vertical profiles proposed in the literature for the region surrounding KKNPP.

Specifically, DA&2010-1 and DA&2010-2 were proposed by [START_REF] Ducellier | Numerical simulation of the mw 6.6 niigata, japan, earthqhake: Reliable input ground motion for engineering purpose[END_REF], in their numerical simulations of the NCOEQ2007 mainshock. Aochi2013 was proposed by [START_REF] Aochi | Finite difference simulations of seismic wave propagation for the 2007 Mw 6.6 Niigata-ken Chuetsu-Oki earthquake: Validity of models and reliable input ground motion in the near field[END_REF], for their finite difference analysis of the NCOEQ2007 mainshock and aftershock sequence. Cirella2008 velocity values were employed by [START_REF] Cirella | Rupture process of the 2007 niigata-ken chuetsu-oki earthquake by non-linear joint inversion of strong motion and gps data[END_REF] to perform a waveform inversion of NCOEQ2007 source mechanism. Finally, the profiles at stations NIG004, NIG016, NIG026 belong to K-NET network and NIGH12 to KiK-Net network [START_REF] Nied | K-net, kik-net[END_REF][START_REF] Aoi | Strong-motion seismograph network operated by NIED: K-NET and KiK-net[END_REF]. Those crustal models are reliable at a very low-frequency range (f< 0.5 Hz).

Other geological models

Figure 3 displays the cross-section SC of a local folded geology striking at N 145 • E and located underneath KKNPP. [START_REF] Watanabe | 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[END_REF] employed it in a 2-D finite element earthquake analysis. The folded structure was constructed from a boring and seismic reflection survey [START_REF] Kobayashi | Geology of the Kashiwazaki District[END_REF]. The SC cross-section spans approximately a region 7.6 km wide and 4.8 km. Seven strata are reported, including the seismic bedrock (see Table 1).

The simulation results rendered by [START_REF] Watanabe | 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[END_REF] showed good agreement with the observed strong motion records at Unit 1, explaining the strong ground motion spatial variability within the nuclear facility (due precisely to the first and second boundaries). Those findings are confirmed by [START_REF] Tokumitsu | 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[END_REF]. [START_REF] Tsuda | Modeling 3D Velocity Structure in the Fault Region of the 2007 Niigataken Chuetu-Oki Earthquake with Folding Structure[END_REF] clarified the role of this folded structure by performing a 3-D finite difference earthquake simulation of the NCOEQ2007 mainshock and aftershock sequence. The 3-D velocity structure they proposed is the combination of a broad model, proposed by Japan Nuclear Energy Safety Organization (JNES, see [START_REF] Kamae | Earthquakes, Tsunamis and Nuclear Risks[END_REF]. Although accurate for the regional wave field, JNES did not include the folding structure. Therefore, [START_REF] Tsuda | Modeling 3D Velocity Structure in the Fault Region of the 2007 Niigataken Chuetu-Oki Earthquake with Folding Structure[END_REF] Vs [ m /s] structure, based on the previous studies of [START_REF] Tokumitsu | 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[END_REF], which consists of 7 cross-sections. This local folding model was interpolated to reconnect to the external broad model. Using the horizontal-vertical (H/V) ratio technique [START_REF] Nakamura | A method for dynamic characteristics estimation of subsurface using microtremor on the ground surface[END_REF], [START_REF] Tsuda | Modeling 3D Velocity Structure in the Fault Region of the 2007 Niigataken Chuetu-Oki Earthquake with Folding Structure[END_REF] matched the synthetic H/V with those of the observed records for some sites on the KKNPP, obtaining a good agreement of the predominant frequencies around 0.3 Hz. The numerical simulation for the NCOEQ2007 showed good agreement with the recorded time-histories at KKNPP, in the range 0-4 Hz.

¡ 20000 ¡ 15000 ¡ 10000 ¡ 5000 0 z [ ] Z-Model R-Model1 ZR-Model D& A2010-1 D& A2010-2 Cirella_2018 NIG004 NIG012 NIG016 Aochi2013 (b)
Finally, Gatti et al. (2018a) validated a similar 3-D geological model by comparing the synthetic records with the recorded shakings from two aftershocks ( JMA 4.4 and JMA 4.2) of the NCOEQ2007 (NIED DMC; [START_REF] Obara | A densely distributed high-sensitivity seismograph network in Japan:Hi-net by National Research Institute for Earth Science and DisasterPrevention[END_REF][START_REF] Okada | Recent progress of seismic observation networks in Japan Hi-net, F-net, K-NET and KiK-net[END_REF], confirming that a simple layered 1-D model (Aochi2013 proposed by [START_REF] Aochi | Finite difference simulations of seismic wave propagation for the 2007 Mw 6.6 Niigata-ken Chuetsu-Oki earthquake: Validity of models and reliable input ground motion in the near field[END_REF])

cannot reproduce the ground motion amplification accurately at Unit 1 and the focalization effect due to the folding structure.

Proposed 3-D geological models

In this work, three geological models are considered, obtained from the interpolation of the available 1-D and 2-D geological information described in the previous subsection, integrated with the geological model of reference, provided by the Geological Survey of Japan for the Niigata area [START_REF] Sekiguchi | 3d subsurface structure model of the niigata sedimentary basin[END_REF]. Available information was alternatively combined in each model, aiming at reproducing the regional and site path effects. Namely, the three models are referred to as:

1. Z-Model: this virtual geology corresponds to the model proposed by Gatti et al. (2018a) and Gatti et al. (2018c) Upper-Teradomari Lower-Teradomari 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6. [START_REF] Tsuda | Modeling 3D Velocity Structure in the Fault Region of the 2007 Niigataken Chuetu-Oki Earthquake with Folding Structure[END_REF]. Layer names are reported in Table 1. Reprinted from Gatti et al. (2018a).

to predict the structural response of the KKNPP Unit 7 reactor building. It blends two geological structures: a 3. ZR-Model: this synthetic 3-D configuration is a combination of the previous two, i.e., the folding structure (Z-model) is included in the regional structure of the Niigata basin (corresponding to the R-model). This hybrid geological structure shares the deepest layers (> 6 km) with the R-Model and the surficial ones with the Z-Model.

The main reason why this geology was considered to coincide with the purposes of [START_REF] Tsuda | Modeling 3D Velocity Structure in the Fault Region of the 2007 Niigataken Chuetu-Oki Earthquake with Folding Structure[END_REF] to include the prediction of both the regional wave field (interacting with the Niigata basin-like geology) and the influence of the local geology on it.

The left and right columns in Figure 4 show four West-East cross-sections of the Z-Model and the R-Model, respectively. In Appendix A, the detailed 3-D cuts of the three geological structures in the Niigata region are reported. It is important to note that the Z-model is relatively regular for depths below 5 km (layered velocity structure proposed by [START_REF] Aochi | Finite difference simulations of seismic wave propagation for the 2007 Mw 6.6 Niigata-ken Chuetsu-Oki earthquake: Validity of models and reliable input ground motion in the near field[END_REF], whereas the R-model shows 3-D sharp basin-like structures reaching a depth of 10 km.

Physics-based numerical simulation

3-D Numerical Model

To assess the uncertainty related to the geological and geophysical models proposed in the previous section, the constructed numerical model covers 64 km × 64 km around the earthquake position (see Figure 1). The model characteristics are summarized in Table 2. All three models exploit the very same hexahedral mesh, neglecting the regional topography. Since different geologies characterize the three models, their accuracy varies from model to model. However, all the models are accurate enough to study the influence of the geology for periods > 0.5 . Earthquake simulations were performed using SEM3D (CEA and CentraleSupélec and IPGP and CNRS, 2017), a high- performance software, implementing the spectral element method (SEM, [START_REF] Komatitsch | The Spectral Element Method: An Efficient Tool to Simulate the Seismic Response of 2D and 3D Geological Structures[END_REF][START_REF] Faccioli | 2D and 3D elastic wave propagation by a pseudo-spectral domain decomposition method[END_REF]. This software is based on RegSEM code [START_REF] Cupillard | RegSEM: a versatile code based on the spectral element method to compute seismic wave propagation at the regional scale[END_REF][START_REF] Festa | The Newmark scheme as velocity-stress time-staggering: an efficient PML implementation for spectral element simulations of elastodynamics[END_REF]. The SEM represents a high-order version of the finite element method (FEM), and it is well known for the relatively easy parallelization algorithm associated (domain decomposition over distributed memory computing cores, using MPI protocol [START_REF] Göddeke | Finite and Spectral Element Methods on Unstructured Grids for Flow and Wave Propagation Problems[END_REF]. Those aspects make the 3-D simulation an attractive tool for seismic hazard analysis [START_REF] Olsen | Simulation of 3D elastic wave propagation in the Salt Lake Basin[END_REF][START_REF] Bradley | On-going challenges in physics-based ground motion prediction and insights from the 2010-2011 Canterbury and 2016 Kaikoura, New Zealand earthquakes[END_REF]. The mesh has higher refinement at shallow layers to preserve the numerical accuracy where slowness values are larger. The maximum theoretical frequency propagated can be estimated as:

= , ⋅ 5 ⋅ Δ (3)
with , being the minimum shear wave velocity in the domain, Δ the corresponding element characteristic size and the number of Gauss-Lobatto-Legendre integration points featuring the mesh element (for the SEM, is usually greater than 5).

Given the intricate geologies at stake and due to the adopted integration rule over 3-D hexahedra, a not-honoring approach is adopted to associate the spatially distributed geological properties to the mesh point locations [START_REF] Casarotti | Cubit and seismic wave propagation based upon the spectralelement method: An Advanced Unstructured Mesher for Complex 3D Geological Media[END_REF]. In other words, the geological discontinuities are not directly meshed, but rather interpolated over the GLL grid. This approach has the advantage of simplifying the inclusion of surfaces of complex geometry, yet inducing some spurious amplification due to spatial interpolation. The influence on the earthquake numerical predictions presented in this study was duly checked and it can be neglected, due to the rather low-frequency range considered.

Finally, Gatti et al. (2018a) verified the limited influence of surface topography on the seismic response at KKNPP, which was not therefore included. Appendix B provides some insights on the computational resources employed for each simulation.

As far as the source is concerned, in order to assess the uncertainty related to the geological model, two earthquake point sources were adopted: (1) EQ1, corresponding to the JMA 4.4 (Japan Meteorological Agency) aftershock of the NCOEQ2007 sequence (employed by other authors (Gatti et al., 2018a)) and (2) EQ2, corresponding to the NCOEQ2007 mainshock hypocenter location. Table 3 summarizes the focal mechanism and location of the two sources. These parameters represent the small aftershocks employed to validate the numerical model by Gatti et al. (2018a) against the observed recordings. Statistical analysis on the source's uncertainty, such as by [START_REF] Shang | Point-source inversion of small and moderate earthquakes from p-wave polarities and p/s amplitude ratios within a hierarchical bayesian framework: Implications for the geysers earthquakes[END_REF], are out of the scope of this paper. For the sake of simplicity, the mainshock source (EQ2) was ; ; represent the strike, rake and dip angles, respectively, estimated by F-NET Centroid Moment Tensor solution [START_REF] Kubo | NIED seismic moment tensor catalogue for regional earthquakes around Japan: quality test and application[END_REF].

represents the source risetime.

Event considered as a point source to be coherent with the idea of this paper of studying the effect of epistemic uncertainty on geological and geophysical information. The full moment tensor components were provided by F-NET. More accurate analysis requires an extended source model, whose rupture path highly affected the mainshock pulse-like seismic response at the site (near-source effects shown by [START_REF] Aoi | Source process of the 2007 Niigata-ken Chuetsu-oki earthquake derived from nearfault strong motion data[END_REF]Gatti et al., 2018a). However, the following sections show the interest in considering a mainshock source due to its Source Time Function (STF) featured by large rise-time compared to EQ1 (see Table 3).

JMA M 0 [Nm] ; ; [ • ] [km] [s] EQ1 (

Synthetic seismic response of the Niigata region

Before performing statistical analyses on the post-processed results (see Section 4), a validation task is presented hereafter. In doing so, synthetic time-histories are compared to recorded signals for EQ1, an aftershock widely recorded by the down hole arrays and surface seismometers located within KKNPP (see Figure 5). The site is instrumented with two recording networks, including four borehole arrays [START_REF] Gatti | Investigation of the earthquake ground motion coherence in heterogeneous non-linear soil deposits[END_REF], with KSH in Figure 5 representing a free-field control point for this analysis (Gatti et al., 2018a).

Given the accuracy of the numerical models, all the synthetic and recorded time-histories were band-pass filtered in

Uncertainty quantification on the synthetic ground motion

In this section, the uncertainty of the numerical earthquake prediction related to the chosen digital geological model is assessed, exploiting EQ2 source. The outcomes of each simulation, obtained by plugging in the three different geological models rendered in Figure 4, are compared in terms of trends and residual values, following the GMPE framework [START_REF] Rodriguez-Marek | A Model for Single-Station Standard Deviation Using Data from Various Tectonic Regions[END_REF]. Compared to GMPE, numerical earthquake simulations yield a more accurate prediction of site-specific seismic response, which is generally difficult to cast into non-linear regressions,

given their intrinsic uniqueness. Moreover, in this paper, up to 261120 virtual stations were deployed across a square grid, allowing a significant statistical analysis in the post-processing phase.

Sensitivity of physics-based simulations to geological structure

For engineering purposes, the pseudo-spectral acceleration (5% damping) is usually adopted as an insightful IM characterizing the local seismic response. Figure 7 shows the contour values of ( = 1 ) for each model (Figure 7a, 7c and 7e), along with the respective contour maps of the thickness of the layer with < 1500 m/s (Figure 7b, 7d, 7f).

The difference between the synthetic seismic response obtained with R-Model (Niigata basin) and the one obtained with the Z-model (layered regional geology) is evident, which leads to infer that this discrepancy relies on the thickness of surficial soil layers. Compared to the Z-Model and the ZR-Model, which share the same horizontal layering close to the surface, the R-model induce a much more heterogeneous distribution.

Z-model seemingly better fits the recorded response (see Figure 6), whereas including the Niigata basin geology into the earthquake simulation (R-model) over-amplifies the response at shorter periods. This statement is supported by Figure 8a, which shows the response spectra for the three geological models considered. It is observed that Z-Model (layered regional geology and local folding structure) provides the lowest values overall. On the other hand, Figure 8b reports the spectral ratio * ∕ * * between two alternative synthetic responses, * and * * , respectively. The spectral ratio assesses the relative amplification -at several natural periods -of two seismic responses predicted by considering two different geological models. The influence of the geological structures interacting with the impinging wave field is highlighted according to their characteristic dimensions (e.g., layer thickness) relative to the incident wavelengths. In the period range, 0.5 s < < 5 s, the ratios range between 0.5 and 1.5, with and exhibiting the far-most comparable response in the range 1.5-3.3 s, where the significant impact of the local folding geology is evident on the synthetic ground motion at KKNPP.

The simulated seismic response was compared to the Abrahamson et al. ( 2016) GMPE, suitable for the seismic context at stake. Figure 9 displays this comparison, for values at = 1 s (all stations in Figure 9a and geometric means in Figure 9c) and = 3 s (all stations in Figure 9b and geometric means in Figure 9d) along the source-to-site distance . ̄ represents the geometric mean response spectra [START_REF] García-Fernández | Modelling Pan-European ground motions for seismic hazard applications[END_REF] expressed as: (d) ).

̄ ( ) = 1 ∑ =1 log ( ) (
with being the response spectra at the station . For each model, the synthetic values (both at = 1 s and = 3 s) are quite dispersed, compared to the ± (i.e. plus/minus one standard deviation) margins indicated by the Abrahamson et al. ( 2016) GMPE, for a reference ,30 = 400 m/s (approximately corresponding to the lowest value in the R-model, see Figure 4). However, values computed at the KKNPP (KSH-SG1 station in Figure 5, referenced by a star in Figure 9a and Figure 9b) for R-model are higher than ̄ , at the same distance from the source. On the contrary, despite the fact that the smallest value in Z-model is 700 m/s (the Nishiyama layer in Table 1) the corresponding value computed at the KKNPP is in very good agreement with the ̄ prediction of the [START_REF] Abrahamson | BC Hydro ground motion prediction equations for subduction earthquakes[END_REF] GMPE. This inconsistency can be explained by the fact that the GMPEs generally lack an explicit description of 3-D site effects, that cannot be effectively condensed into the ,30 value.

It is interesting to notice that trends in Figure 9c and Figure 9d can be adopted as site-specific and scenario-specific GMPE predictions. A stable decaying trend is observed within the 0-35 km distance range, where the number of stations is statistically representative. The thinner solid lines in Figure 9c and Figure 9d same way along with the source-to-site distance, but they consistently drop below 0.4 at longer period (see Figure 10b). A possible explanation comes from the impact of the deep geological structure of the Niigata basin (R-model) on the long-period ground motion prediction ( ( = 3 s)), which is remarkably different from the one estimated with the layered configuration in the Z-model. However, the influence of the complex yet deep geological structure seemingly vanishes at shorter periods (e.g., at = 1 s), where values are mainly affected by the folding structure located within the first 5 km below the surface.

Hereafter, the synthetic seismic response estimated for the three geological models is analyzed in order to estimate the residual of the synthetic strong ground motion prediction.

Estimation of within-event variability

Several studies (e.g., [START_REF] Atkinson | An empirical perspective on uncertainty in earthquake ground motion prediction[END_REF][START_REF] García-Fernández | Modelling Pan-European ground motions for seismic hazard applications[END_REF] estimated the standard deviation of the withinevent residual at a certain source-to-site distance as: with the source-to-site distance. This constant trend is insignificant for = 3 s. Since those models share the same shallow horizontal layering (including the folding), is mainly related to shallow geology, indicating once again its major impact compared to deep geology. Exploiting the rare advantage of having many stations close to the source for a single event, Figure 11 proves that is approximately independent of the distance, even in the case of R-Model where there is a complex layering (non-horizontal). In many cases, the value of is considered to be independent on the distance to the source [START_REF] Atkinson | An empirical perspective on uncertainty in earthquake ground motion prediction[END_REF]. However, recent studies proposed a functional dependence of on the distance from the source (see [START_REF] Boore | NGA-West2 equations for predicting PGA, PGV, and 5% damped PSA for shallow crustal earthquakes[END_REF], although this variation occurs at distances larger than 80 km.

( ) = √ √ √ √ 1 -1 ∑ =1 log ( ) ̄ ( ) 2 (5)

Influence of the site effects on the within-event variability

To explain the variability of the within-event residuals (

), its sensitivity with the respect to the local site conditions is hereafter analyzed. In doing so, is related to several parameters such as: (1) the harmonic-mean shear velocity in the first m below the surface ( , ) and ( 2) the thickness of the equivalent basin-like structure ( ). Those parameters are used in the recipes proposed by Eurocode 8 (EC8, CEN, 2004) and they are currently being discussed in the framework of the EC8 site categorization update [START_REF] Paolucci | Checking the site categorization criteria and amplification factors of the 2021 draft of Eurocode 8 Part 1âĂŞ1[END_REF]. , is often computed for = 30 m and = 100 m. The thickness of the equivalent basin-like structure is instead defined based on the shear-wave velocity value of the engineering bedrock, such as 800 , corresponding to the depth at which the engineering bedrock is found, i.e. the geological layer at which the shear velocity ranges around 800 m/s (according to [START_REF] Ansal | Source and site factors in microzonation[END_REF].

In this study, 1500 is also considered, in order to assess the role of deeper geological layers (such as the Shiiya and

Upper-Teradomari layers in the folding structure, see Figure 3). Table 4 shows the correlation between the R-Model parameters for each mesh point location at surface. As expected, ,30 and ,100 are well-correlated with = 0.98.

Similarly, the thicknesses 800 and 1500 are correlated with = 0.78. ,30 and ,100 and 800 and 1500 are poorly correlated and they can be used as complementary parameters to predict better and reduce . Figure 12 shows values as a function of the four parameters ,30 , ,100 , 800 and 1500 respectively. Positive indicate that the synthetic ground motion at the station of interest is larger than the geometric mean prediction. Negative values, on the contrary, point out that the synthetic IM at the station underestimate the geometric mean prediction.

Adopting the R-model, the estimation (for both = 1 s and = 3 s) is positive but it rapidly decays with ,30

and ,100 values below 800 m/s, whereas it remains approximately constant and close to -1 for larger shear-wave velocity values. Numerical simulations seemingly overestimate the ̄ values across softer sediments, highlighting the presence of 3-D site effects that the GMPEs can barely picture. Moreover, the trend in Figure 12 are in agreement with the EC8 classification and with the GMPE guidelines (e.g., [START_REF] Abrahamson | Attenuation of long period strong ground motions[END_REF][START_REF] Lee | Empirical equations describing attenuation of peaks of strong ground motion, in terms of magnitude, distance, path effects and site conditions[END_REF] that usually classify a site according to its ,30 , being a soil deposit for ,30 < 800 m/s and engineering bedrock for ,30 ≥ 800 m/s. In addition, this result supports the choice of 800 to condense the basin-like site effects.

On the contrary, no clear trend is observed for the estimated at = 1 s neither with respect to 800 nor with the respect to 1500 (see Figure 12). At longer period (i.e., = 3 s), as suggested by decaying trend approaching small 800 or 1500 , the residual diminishes for lower 1500 , highlighting the influence of harder and deeper rock layers on the site effect prediction.

Figure 13 shows the absolute correlation coefficient | | between and each parameter ,30 , ,100 , 800 and 1500 , at = 1 s (blue histogram) and = 3 s (orange histogram). In this figure, the same correlation is computed for the Z-model for the sake of comparison. | | with respect to 800 and 1500 almost doubles at long period ( = 3 s) highlighting the importance of site effects on the residual estimation. In the case of Z-Model (Figure 13b), there is no correlation between and ,30 and poor correlation with ,100 , due to the horizontal layering and the coarse geological model. Both parameters ,30 and ,100 are uninformative at this scale and resolution. However, 800

and 1500 represent complementary site-characterization parameters to ,30 . Finally, the intensity of the ground motion is affected by the pattern of the source. Figure 7 showed that spatial variability at surface depends on the station back-azimuth. Figure 14 For the ZR-model, featured by a horizontally layered regional geology (see Figure 14a), the clustering based on the back-azimuth highlights a maximum difference of up to 40% between clusters. When considering the R-Model (see Figure 14b), those discrepancies between cluster are less evident. For back-azimuths between 7 • and 127 • , the high values are seemingly related to the complex geology in this zone. This is congruent with several studies (such as [START_REF] Ripperger | Variability of Near-Field Ground Motion from Dynamic Earthquake Rupture Simulations[END_REF][START_REF] Vyas | Distance and Azimuthal Dependence of GroundâĂŘMotion Variability for Unilateral StrikeâĂŘSlip Ruptures[END_REF] that also proved the variation of as a function of the source pattern effects.

In this study, point sources are considered, although the back-azimuth dependency could be higher for extended-fault scenarios due to directivity effects.

Estimation of between-event variability

As far as is concerned, the extensive number of monitoring points allows to draw solid statistics in a large source-tosite distance range. In this study, the quota of epistemic uncertainty linked to the source mechanism is not investigated. Therefore, can be indirectly estimated by disentangling it from the within-event variability. Figure 15 shows the estimated for both ( = 1 s) and ( = 3 s). However, the fact that the three cases considered herein display the same source pattern and location affects its estimation [START_REF] Atik | The Variability of Ground-Motion Prediction Models and Its Components[END_REF]. Concerning the variability of the betweenevent residual , one easily observes that the variance of the total residual 2 is the sum of the variances of betweenevent and within-event residuals, which are modeled as independent random variables. Therefore, the between-event variance is estimated as:

2 = 2 -2 (6)
with 2 estimated in Section 4.2 from the synthetic time-histories. Therefore, the first step is to compute the modelto-model deviation, inferred directly from the * ∕ * * ratio, according to the following expression:

Δ * ∕ * * = Δ * -Δ * * (7)
where Δ , for any model and period, is considered as distributed as centered log-normal random variable. With this assumption, the variance of Δ * ∕ * * , namely 2 * ∕ * * is computed as: is very low. This aspect, along with the high correlation between both models (see Figure 10), indicates that, even if they do not share similar geologies at depth, the two models are almost interchangeable for this very earthquake.

2 * ∕ * * = 2 * + 2 * * -2 ⋅ * , * * ⋅ * ⋅ * * (
The second step to estimate the variability of consists into assembling a system of the three equations in the form of Equation ( 8), considering , and as unknowns. The system is iteratively solved by exploiting the outcome of the three earthquake simulations. Finally, is computed by exploiting Equation (6). Figure 17 depicts the estimation along the source-to-site distance (not taking into account the variability that could come from the model or the parameterization of the source event). Figure 17 shows that the three curves are slowly varying along the source-to-site distance. Despite some minor fluctuation between 20 and 24 km, is practically constant for all the three models at = 3 s. Moreover, the curves referring to the ZR-Model and the Z-Model are very similar for both = 1 s (slowly decaying values with increasing distance from the source) and = 3 s (lower deviation, constant along the source-to-site distance). The deviation associated to the R-model is almost insensitive to the natural period ( ≈ = 0.45) except for its drop to lower values at short distances and long period, arguably due to the rather constant geological properties (compared to the layered medium in Z-model) at the surface depth.

Conclusions

This study compares the synthetic outcomes of three physics-based numerical models, each featuring a different yet plausible geological structure for the region at stake. A realistic test case (the 2007 Niigata earthquake) is studied, so to prove the sensitivity of high-fidelity earthquake simulators to input geological information. The comparison proved that complex geological features (the Niigata regional basin plugged into the R-Model), compared to horizontally layered deep geology and folding geology close to the surface (Z-model and ZR-Model) might induce less than 10% difference on the surface seismic response (in terms of values), leading to the possibility of simplifying the geological medium to reduce the computational costs). Moreover, physics-based simulations were employed to study the geology's influence on the spatial variability of the ground shaking.

Considering the 261120 virtual stations covering the region of interest, the paper adopted the GMPEs formulation in Equation ( 1) to infer synthetic site-specific ground motion prediction variability. From the analysis, we evinced that (i.e. the within-event residual standard deviation) mostly depends on shallow layers properties, increasing with the geological complexity. Furthermore, the two models sharing the surficial site conditions but different ones at deep layers (Z-Model and ZR-Model) show a remarkable difference in the estimated between-event residual by an almost constant factor (with source-to-site distance). This suggests that records from different tectonic regions but similar site conditions (basin-like structures) can be integrated into a site-specific GMPE, by applying a correction factor that considers the diverse deep geology.

The impact of complex layering on (i.e. the between-event residual standard deviation) is reflected in the R-model results, where higher is obtained than in simpler geology models. Shallow geology is proven as a significant factor in determining the average trend and variability of the synthetic seismic site response. At long period ( = 3 s), the between-event variability is seemingly highly affected by the choice of the geological model, in a source-to-site distance range of 0-25 km.

The outcome of the numerical exercise suggests the possibility of disentangling the contribution of shallow and deep geological structures on the overall synthetic seismic prediction at the site.

Finally, we unveiled a relation between the within-event residuals ( ) and the properties of the sedimentary layers, such as basin thickness and harmonic shear-wave velocity. Those two parameters are seemingly equally correlated to the residual . In particular, the depth of the engineering bedrock is a rather easily accessible information, inferred from regional geological maps or obtained by non-invasive geophysical tests. Therefore, either 800 or 1500 are suitable to reduce the variability of the GMPE prediction, especially at long period. In this study, the influence of the source pattern and the back-azimuth on both and variability is assessed.

For the present case study, the results of this numerical analysis are congruent with several previous studies on real records. This shows that the numerical simulation from a well-calibrated model can be integrated with real records to create synthetic GMPEs, or complete the database of a region to produce more accurate GMPE. The previous results

were found with a frequency band 0-2.1 Hz and not including the topographical effects. It is possible that at higher frequencies, the local conditions are more relevant, especially in the value of . However, the results found in this simulation are still valid.

A. Details of the geological models

B. Performance of the numerical analysis

Figure 1 :

 1 Figure 1: Digital Elevation Model of the region considered in the analysis, including the Niigata basin and the surroundings of the Kashiwazaki-Kariwa Nuclear Power Plant (KKNPP). Yellow dashed horizontal lines represent the traces of the vertical clips shown in Figure 4.

Figure 2 :

 2 Figure 2: 1-D velocity models of the different boreholes from the region around KKNPP. The triangles are located at the maximum borehole depths. a) Compressional wave speed profile; b) Shear wave speed profile.

Figure 3 :

 3 Figure 3: Cross-section (SC) of the folded geology below KKNPP (Unit 1 and Unit 5) proposed by Watanabe et al. (2009); Tsuda et al. (2011). Layer names are reported in Table1. Reprinted fromGatti et al. (2018a).

  regional sub-horizontally layered geology inspired by[START_REF] Aochi | Finite difference simulations of seismic wave propagation for the 2007 Mw 6.6 Niigata-ken Chuetsu-Oki earthquake: Validity of models and reliable input ground motion in the near field[END_REF] (labeled as Aochi2013 in Figure2) and a model obtained from the extrusion of SC cross-section in Figure3. The two geological structures are joined by linear interpolation. The Aochi2013 model was calibrated to reproduce the regional wave field in the 0-0.5 Hz frequency range. The hybrid model, including the folding, was instead verified in the 0-5 Hz frequency range, for two NCOEQ2007 aftershocks, at several locations within the KKNPP site (the borehole arrays deployed within the nuclear site) and at several KiK-Net and K-NET stations.2. R-Model: this model was proposed by the Geological Survey of Japan (GSJ,[START_REF] Sekiguchi | 3d subsurface structure model of the niigata sedimentary basin[END_REF]. It is an improvement of the previous NIED model (National Research Institute for Earth Science and Disaster Prevention) proposed by[START_REF] Fujiwara | A subsurface structure modeling of whole of Japan for strong-motion evaluation[END_REF] for the Niigata area, based on several seismological observations. The depth of layer boundaries is provided on a grid mesh of approximately 0.5 × 0.5 km. The model contains 50 layers above the Moho (Mohorovicic discontinuity).[START_REF] Aochi | Finite difference simulations of seismic wave propagation for the 2007 Mw 6.6 Niigata-ken Chuetsu-Oki earthquake: Validity of models and reliable input ground motion in the near field[END_REF] tested and compared this geological structure with two other 3-D models of the Niigata basin and reproduced the NCOEQ2007 mainshock time-histories at several KiK-Net (NIGH11) and K-NET stations (NIG019 and NIG017) in the frequency range 0.1-0.5 Hz. The authors specified that the synthetics at near-field soil sites were quite realistic, but the model needed improvements and a finer shallow structure.

Figure 4 :

 4 Figure 4: West-East cross-sections at different latitudes (reference system: Tokyo CS VIII EPSG:30168 system), showing the contours. Figure 1 shows the position of the cross-sections in the map. Left column belong to the Z-Model and right column to the R-Model.

Figure 5 :

 5 Figure 5: Map of the KKNPP site. The locations of downhole arrays of seismometers are reported: KSH represents the Service Hall array (four devices, SG1-SG4, reaching depth of -250 m G.L.), which entirely recorded the NCOEQ2007.

Figure 6 :

 6 Figure 6: Average time-frequency Envelope GoF (EG) and Phase GoF (PG) (proposed by Kristekova et al., 2009) for recording stations at KKNPP, during EQ1 (aftershock). Synthetic and recorded time-histories were band-pass filtered in the 0.1-2.8 Hz frequency range. GoFs were computed between each numerical model (R-model (blue), Z-model (orange) and ZR-model (green) and recorded signals. The histograms collecting EG and PG GoFs are reported on each axis.

Figure 7 :

 7 Figure 7: Contour plot of values at = 1 s (5% critical damping) for the three different models: (a) Z-model, (c) R-model, (e) ZR-model. Contour maps of the thickness of the geological layers with < 1500 m/s for the three models: (b) Z-model, (d) R-model, (f) ZR-model.

Figure 8 :

 8 Figure 8: Comparison of the response spectra at the KKNPP location (KSH -37.4 • N, 138.6 • E) for the three geological models considered. (a) measures the geometric average of the horizontal components. (b) measures the ratio * ∕ * * between two alternative synthetic response * and * * , corresponding to different geological models. Orange area covers the frequency band highly affected by the applied filter ( < 1∕).

  Figure 9: values obtained at surface for the three physics-based simulations at: (a) and (c) = 1 s; (b) and (d) = 3 s. (a-b) show the synthetic response at each station, compared with the GMPE proposed by Abrahamson et al. (2016) (computed for a reference ,30 = 400 m/s. The geometric mean ̄ , for the three physics-based simulations, is depicted with a solid thick line in (c-d), whereas the respective ± margins are marked by thinner dashed lines. Stars indicate the free-field response at KSH station, located within the KKNPP.

Figure 10 :

 10 Figure 10: Correlation coefficient , of the values, computed between the th and th numerical models as a function of the source-to-site distance . (a) Correlation coefficient of ( = 1 s); (b) correlation coefficient of ( = 3 s).

Figure 11 Figure 11 :

 1111 Figure 11 displays the estimation of for the present study, at different source-to-site distances and for = 1 s and = 3 s. For Z-model and ZR-model and for both natural periods, Figure 11 depicts a rather constant trend along

Figure 12 :Figure 13 :

 1213 Figure 12: Within-event residual computed for the R-model, as function of ,30 , ,100 , 800 and 1500 . Top panels: computed for ( = 1 s). Bottom panels: computed for ( = 3 s).

Figure 14 :

 14 Figure 14: Standard deviation of the within-event residual as a function of the source-to-site distance . values at each station are clustered based on the station back-azimuth. (a) ZR-Model; (b) R-Model.

Figure 15 :Figure 16 :

 1516 Figure 15: Between-event residuals for (a) ( = 1 s) and (b) ( = 3 s) as functions of the source-to-site distance .

Figure 17 :

 17 Figure 17: Standard deviation of the between-event residual computed along source-to-site distance , for = 1 s (a) and = 3 s (b).

Figure 18 :Figure 19 :Figure 20 :

 181920 Figure 18: View-clips of the Z-Model. A) Superior view at = -1500 m. B) Longitudinal cut by the South-North direction (y-Axis). C) different slips by the East-West (x-axis) direction.

  This work was granted access to the HPC resources of CINES under the allocation 2018-A0040410444, 2019-A0060410444, 2020-A0070411083 and 2020-A0080410444, made by GENCI (Grand équipement national de calcul intensif). Computations were also performed using HPC resources allocated by the Mésocentre Moulon, the super-computer of Cen-traleSupélec and École Normale Supérieure de Paris Saclay, Paris Saclay University.Given the scalability properties of SEM3D (https://www.researchgate.net/publication/349254101_SEM3D-High-res seismic_wave_propagation_modelling_from_the_fault_to_the_structure_for_realistic_earthquake_ scenarios_GENCI_Allocation_A0080410444/stats), each run (30 s earthquake simulation, on the template model) took approximately 80 minutes (≈ 960 hours CPU-time) on 720 MPI cores Intel Xeon Gold 6230 20C @2.1 GHz Cascade Lake.C. Numerical model validationIn Figure21, the synthetic and recorded accelerograms and Kristekova's GoF at KSH (SG4) are compared, considering each numerical model (R-model, Z-model and ZR-model respectively). As expected, the synthetic ground motion rendered at the KKNPP differs remarkably from model to model. Moreover, Z-Model displays a shorter coda-wave compared to the R-Model, where the scattering effect caused by the more intricate geological structure increases the effective propagation path and rendering large-duration time-histories.

Figure 21 :

 21 Figure 21: Synthetic acceleration time-histories (EW component, filtered between 0-2.1 Hz) at KSH location (device SG4, G.L. -250 m) for the three models for EQ1. Red time-histories refer to earthquake recording (TEPCO, 2007), whereas black time-histories refer to the outcome of the numerical simulation. (a) R-model; (b) Z-model; (c) ZR-model.

Table 1

 1 Geological properties of the folding structure underneath KKNPP. and are the compressional-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 another.

	Layer	[m/s]	[m/s]	[kg/m 3 ]
	Nishiyama	700	1900	1700
	Shiiya	1200	2200	2100
	Upper Teradomari	1700	3300	2300
	Lower Teradomari	2000	4200	2400
	Nanaya	2000	4600	2500
	⋆⋆ Green tuff	2600	5200	2600
	Seismic bedrock 2600->2750	5200	2600

Table 2

 2 Parameters of the template numerical model employed in this paper. and represent the number of hexahedral elements and degrees of freedom (DOF) respectively. Δ is the minimal element characteristic size, is the number of Gauss-Lobatto-Legendre integration points per element size, is the maximum frequency that the model can theoretically propagate.

	Model	[1]	[1]	Δ	[m]	[1]	[m/s]	[Hz]
	Z	≈ 4.02⋅ 10 6 ≈ 1.50 ⋅ 10 9	139	5 ×5×5	700	5
	R/ZR	-	-	-		-	400	2.1

Table 3

 3 Summary of the aftershock parameters employed in this analysis.

Table 4

 4 Correlation matrix of the four analyzed site parameters for the R-model: ,30 , ,100 , 800 and 1500 .

		,30	,100	800	1500
	,30	1.0	0.98	0.40	0.40
	,100		1.00	0.50	0.46
	800			1.00	0.78
	1500	sym			1.00

  confirms this assumption, displaying the value of by clustering each station based on its back-azimuth.

			ZR-Model						R-Model
	0.8						0.8
					187 o t o 247 o						187 o t o 247 o
	0.7				247 o t o 307 o	0.7					247 o t o 307 o
					307 o t o 7 o						307 o t o 7 o
	0.6				7 o t o 67 o	0.6					7 o t o 67 o
					67 o t o 127 o						67 o t o 127 o
	0.5				127 o t o 187 o	0.5					127 o t o 187 o
					All							All
	0.4						0.4
	0.3						0.3
	0.2						0.2
	0.1						0.1
	0.0						0.0
	18	20	22	24	26	28	30	18	20	22	24	26	28	30
			(a)							(b)

Acknowledgements

The research reported in this paper has been supported in part by the SEISM Paris Saclay Research Institute. 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, from the FY 2018 Postdoctoral Fellowship issued by the Japanese Society for Promotion of Science (JSPS) International Fellowships for Research in Japan, hosted by Prof. H. Kawase at Disaster Prevention Research Institute (DPRI) and by the digital geological data provided by Prof. H. Sekiguchi, as well as from valuable strong motion records during the 2007 Niigataken Chuetsu-Oki earthquake disseminated by TEPCO (Tokyo Electric Power Company).