Regional physics-based simulation of ground motion within the Rhȏne Valley, France, during the MW 4.9 2019 Le Teil earthquake

In this paper we introduce the 3D physics-based numerical simulations (PBS) of ground motion during the Nov 11, 2019, MW 4.9 Le Teil earthquake, which occurred in a low-to-moderate seismicity area in the South-East of France, within the Rhône river valley, which hosts several operating nuclear installations. The numerical code SPEED, developed at Politecnico di Milano, Italy, was used to produce the PBS. After introducing the criteria to construct the numerical model, based on the relatively limited data available, a numerical convergence test was made to identify the frequency range for accurate simulations. Furthermore, the performance of the numerical results against the available strong motion records was assessed quantitatively using Goodness-of-fit (GoF) measures. According to the GoF scores, a good-to-excellent agreement was found on the horizontal components up to 8 Hz, showing that, even without a very detailed 3D numerical model of the medium, that would imply detailed knowledge of the basin shape, of the bedrock-to-basin impedance ratio, as well as of the damping ratio in the basin and its dependence on frequency, the PBS may provide realistic broadband predictions of earthquake ground motion. Nevertheless, as shown by the poorer performance on the vertical component, the high-frequency limitations of PBS, also in relation to the energy radiated by the kinematic source model, is still an issue to be carefully addressed. In spite of these limitations, the results obtained in this work demonstrate that PBS, if suitably calibrated and validated, can be either an alternative or a useful complement to empirical ground motion models, especially in those cases where the region- and site-specific features of ground shaking, including near-source conditions, are typically not accounted for by ergodic empirical models, such as for the seismic risk evaluation of large urban areas and/or of strategic structures, infrastructures and industrial plants.


Introduction
On November 11, 2019 a M W 4.9 earthquake occurred in a densely populated low-to-moderate seismicity region of South-Eastern France, within the lower Rhône Valley, very close to the town of Le Teil, at relatively short distance from two nuclear power plants (NPP), i.e., Cruas (at an epicentral distance R epi of 15 km) and Tricastin (R epi = 24 km). Owing to the proximity of the NPPs to the epicenter of the earthquake, the Le Teil earthquake (see Fig. 1) has raised further public attention that seismic hazard may be a relevant contributor to the overall risk of NPPs, as dramatically documented by the accidents at Kashiwazaki-Kariwa (during the July 16, 2007, Chuetsu earthquake, M W 6.6) and Fukushima (during the March 11, 2011, Tohoku earthquake, M W 9) NPPs in Japan.
According to IAEA safety standards (IAEA-SSG9 2022) both empirical and direct simulation methods can be used to estimate vibratory ground motions, within either probabilistic or deterministic seismic hazard assessment. Empirical Ground Motion Models (GMMs) represent the standard approach for ground motion characterization in a probabilistic framework, suited to account for different sources of aleatory and Fig. 1 Seismotectonic map of the region hit by the November 11, 2019 M W 4.9 Le Teil earthquake. The orange and purple circles are instrumental and historical earthquakes, for which moment magnitude (M W ) and maximum macroseismic intensities (I) are given, respectively; the yellow hexagons denote the Nuclear Power Plants (NPP) in the region. The brown lines are potentially active faults mapped in the BDFA catalogue (https:// bdfa. irsn. fr/), while the La Rouvière Fault (LRF), activated by the Le Teil earthquake, is in red. Blue triangles are the stations of the RESIF and EDF networks epistemic uncertainties, but they are also used as a standard for deterministic seismic hazard assessment. However, because of their ergodic nature, classical GMMs do not provide quantitative estimates of the region-and site-specific features of earthquake ground motion, unless empirical non-ergodic adjustments are considered (e.g., Biro and Renault 2012;Ameri et al. 2017) or fully non-ergodic models are implemented in the considered region (e.g. Landwehr et al. 2016;Sung et al. 2022).
As an alternative approach to GMMs, 3D physics-based numerical simulations (PBS) of seismic wave propagation account for the seismic source, the propagation path and the amplification effects related to the site-specific shallow geology (e.g., Paolucci et al. 2018;McCallen et al. 2021). They are becoming more and more appealing as the performance of computer codes is growing exponentially and their use is particularly appropriate in case of complex geological configurations, coupled with near-source conditions, cases that are poorly constrained in GMMs due to small amount of recordings. Quoting IAEA-SSG9, the PBS procedures "might be especially effective in cases where nearby faults contribute significantly to the vibratory ground motion hazard at the site and/or where the existing empirical data are limited (e.g. on the hanging wall of a nearby fault)".
Validation of 3D PBS against recorded weak or strong ground motions is one of the key propaedeutic activities to ensure that the different input elements of PBS, namely, the seismic source and the 3D velocity model, are suitable to reproduce the recorded motions, at least up to a prescribed frequency limit.
In the framework of the SIGMA-2 Project (https:// www. sigma-2. net/), funded by different industrial partners that operate in the nuclear energy sector, a benchmark on different simulation approaches for earthquake ground motion prediction was organized, with reference to the Le Teil earthquake (El Haber et al. 2022). The primary goal of the benchmark was to validate and explore the potential of different ground motion simulation techniques in predicting ground motion in a low-to-moderate seismicity area, where the description of the seismic wave propagation medium is limited, the fault geometry and activity are poorly known and the earthquake records are rare.
In this paper, the simulations are carried out using the spectral element code SPEED (http:// speed. mox. polimi. it/, Mazzieri et al. 2013), which has been extensively used in the recent past to perform PBS validated on different real earthquakes (Paolucci et al. 2015;Evangelista et al. 2017;Infantino et al. 2020;Sangaraju et al. 2021) and to construct a prototype dataset of simulated broadband accelerograms with the aim of complementing recordings datasets, still relatively sparse in such near-source conditions .
After a brief overview of the case study in Sect. 2, the 3D numerical model is introduced in Sect. 3, while the verification and numerical convergence tests are discussed in Sect. 4. In Sect. 5, simulated ground motion is successfully compared with the recorded one on a broad frequency range, up to about 8 Hz, where the convergence tests have shown that the accuracy of numerical wave propagation is preserved. Furthermore, Goodness-of-Fit tests show good to excellent scores. Since the main role of 3D PBS is to highlight region-and site-specific features of earthquake ground motion that cannot be resolved by the ergodic empirical GMMs, and that may lead to biased estimates for seismic hazard assessment, in Sect. 6 the main findings related to the 3D site amplification features in the Rhône Valley are summarized and compared with 1D approaches for site response estimation.

Case study: the Le Teil M W 4.9 earthquake
In spite of its moderate magnitude, the earthquake caused different degrees of damages to approximately 900 residential and several public buildings in the municipality of Le Teil, located at 4 km from the epicenter, going from light cracks in the walls to total collapse (Ritz et al. 2020). A maximum macroseismic intensity (EMS98 scale, Grünthal et al. 1998) I max = VII-VIII was estimated for Le Teil municipality (Schlupp et al. 2021;Sira et al. 2020). The economic losses induced by the earthquake have been estimated at around 200 MEUR for private property and at around 12 MEUR for communal properties (AFPS 2021).
The region is characterized by low-to-moderate seismicity, with instrumental earthquakes of maximum magnitude ranging between 3 and 4 (see orange circles in Fig. 1, from SI-HEX (Cara et al. 2015) updated catalogue, https:// www. franc eseis me. fr/). The most significant historical earthquakes in the region (as indicated by purple dots in Fig. 1) occurred south of Le Teil in 1773, 1873 the 1923, with I max up to VII MSK (SISFRANCE database, www. sisfr ance. net). The August 8, 1873 earthquake, at around 8 km southward from Le Teil, was the largest shock ever felt in this region, with an estimated M W of around 4.1 and a focal depth of about 3 km (FCAT catalogue, Manchuel et al. 2018). An earthquake was located near Le Teil in November 1923, with an inferred M W of around 3 and I max = IV MSK.
From a seismotectonic point of view, the epicenter of the Le Teil earthquake is located at the boundary between the Massif Central crystalline basement and the sedimentary basin of South-Eastern France bordering the Alps mountain range. The tectonic evolution of this region was marked by several deformation phases since 200 million years (Ma), which have produced a complex structural pattern in a compressional stress regime with around 100-km-long system of faults (i.e., the Cevennes Fault System-CFS) striking NE-SW and dipping to the southeast (Delouis et al. 2021;Ritz et al. 2020).
The earthquake was generated by the seismic rupture of a segment of the La Rouvière fault (LRF, see red line in Fig. 1), which is located at the North-Eastern part of the CFS. The LRF was not identified as a potentially active fault in the Database of Potentially Active Faults for Metropolitan France-BDFA (https:// bdfa. irsn. fr/, Jomard et al. 2017), but it was already listed on the geological map of the Aubenas area (Elmi et al. 1996). The 8 km-long La Rouvière fault is oriented NE-SW (azimuth from N030 to N050), it dips steeply to South-East and is located between, and parallel to, the Saint Remèze fault (part of the Cévennes fault) to the North-West and the Marsanne fault to the South-East. The latter two, contrarily to the LFR, were identified as potentially active faults in the BDFA.
Geodetic, seismological and field data indicate a rupture area of about 4 km × 1.7 km, characterized by a reverse focal mechanism, with a hypocenter (44.521°N; 4.669°E) located at just 1 km depth from the ground surface Causse et al. 2021). Ritz et al. (2020) show evidence of surface fault rupture with a permanent uplift up to 15 cm on the fault hanging wall, which is rather uncommon considering the magnitude and for this region. Such a shallow focal depth is unusual for an earthquake of tectonic origin and it is typically associated with earthquakes of anthropogenic nature. Referring to the wellknown cases of seismic events in Groningen (the Netherlands) and in Oklahoma (USA) induced by gas extraction and fluid injection activities, respectively, the analysis of the recorded datasets (publicly available at CESMD, https:// www. stron gmoti oncen ter. org/, and KNMI, http:// rdsa. knmi. nl/ datap ortal/, portals) indicates that most induced earthquakes, with magnitude in the range 2-5, have a focal depth in the range 3-5 km. Based on both in-field observations and numerical simulations, Causse et al. (2021) showed that, although the average source properties of the Le Teil earthquake (stress drop, slip distribution and rupture velocity) were consistent with common deeper earthquakes, the unusually shallow rupture produced exceptional levels of ground shaking in the immediate vicinity of the causative fault. Azimuthal and frequency dependencies of ground motion decay with distance are the object of current research works. The shallow hypocenter, together with the presence of a large limestone quarry located on the hanging wall of the LRF, motivated studies on the causal relationship between the extraction activities and the triggering of the Le Teil earthquake which are still under debate (De Novellis et al. 2020;Liang and Ampuero 2020).
Recordings of the Le Teil event and aftershocks are available from the seismological stations of the RESIF network (Réseau Sismologique et géodésique Français-RESIF http:// seism ology. resif. fr/, 1995) and from the closest stations of the EDF (Electricité de France, the French NPP operator) network. As shown by the blue triangles of Fig. 1, only four stations fall within the area covered by the 3D numerical model (details of these stations are given in Table 1). These stations are the reference with respect to which the simulation results will be tested and validated. Due to the limited number of records at short epicentral distances, PBS can be effectively employed to gain insights into the main features of seismic shaking in the region.

Set-up of the 3D velocity model of the Rhône River Valley
The construction of the 3D subsoil model of the region implied some preliminary analyses to properly identify the main features of the Rhône Valley and of the seismic wave velocity model for both the crustal layers and the sedimentary materials within the valley. In particular, the limited geophysical and geological information at large-scale required to develop a numerical algorithm to shape a preliminary 3D model of the Rhône Valley, constrained on the sparse data made available. Namely, the 3D model of the Rhône Valley shown in Fig. 2 was constructed from numerical processing of the information included in the DEM (Digital Elevation Model) of the area, available at https:// downl oad. gebco. net/, with a resolution of 300 m, further constrained by the sediment depth from available geological cross-sections and by the surface contour of outcropping sediments, as provided by Gélis et al. (2022). For this purpose, an ad hoc algorithm was developed, providing an estimate of the local depth of the Rhône Valley sediments, based on the equilibrium of an elastic and homogeneous membrane fixed at the valley boundaries (i.e., Poisson equation), subjected to a distributed loading inversely proportional to the distance of the point from the boundary. Further details of the procedure are provided in El Haber et al. (2022).
As noticeable from Fig. 2, the valley shape and depth change considerably, being narrow and shallow in the North, close to Cruas NPP, and large and relatively deep in the South, close to Tricastin NPP. The maximum sediment thickness reaches about 700 m near the OGLP station.
The velocity model of the deep crustal layers implemented in the numerical model (Table 2) was borrowed from Causse et al. (2021), who performed a set of numerical simulations of the Le Teil earthquake and characterized for that purpose the 1D structure of the earth crust using seismic noise recorded at temporary seismological stations installed after  Ritz et al. (2020) the earthquake in the fault vicinity. These profiles, in the epicentral area, exhibit soil materials with increasing stiffness from the surface to 1.2 km depth, overlaying less competent deposits (see the geological cross section in Fig. 2). As remarked by Causse et al. (2021), this inversion in the velocity profile is consistent with the geological settings of the area (Elmi et al. 1996) and with information from deep boreholes in the region.
Concerning the sediments, a seismic velocity model was calibrated based on the measured profiles available at the OGLP station (RAP-ID project, Régnier et al. 2010) and at the Tricastin NPP (from local investigations carried out when the NPP was under construction, EDF communication). Other profiles north of Montélimar, around the NPP of Cruas (at CRU1 station), at the border of the basin, were used for verification. Based on this information, a parabolic V S /V P profile was defined as a function of the depth from the topographic surface (z), as follows: For soil density, a constant value ρ = 1.95 t/m 3 was chosen, in agreement with available data. The adopted V S model (Eq. 1) is shown in Fig. 3 (black line), together with the measured profiles and the first layer of the crustal model of Table 2 (dashed brown line). Note that, for sake of simplicity, the velocity profile is homogeneous along horizontal plans in the basin. At the generic point in the basin, the V S profile consists of Eq. (1) until the depth of the basin is reached and then, beyond that depth, the crustal model applies.
Concerning anelastic attenuation properties, for all soil layers, a frequency-dependent quality factor (Q = Q 0 *f/f 0 ) was adopted, with Q 0 = V S /10 and a reference frequency f 0 = 1 Hz.

Kinematic seismic source model
A kinematic representation of the fault rupture process was adopted to model the seismic source of the Le Teil earthquake. In spite of the moderate magnitude, a finitefault modelling was preferred to point-source to provide more realistic ground motion predictions in the near-source region.
Among the studies devoted to the inversion of a kinematic model of the seismic source (Delouis et al. 2021;Cornou et al. 2021;Ritz et al. 2020;De Novellis et al. 2020;Mordret et al. 2020), the one proposed by Cornou et al. (2021) was adopted, in agreement with the partners of the SIGMA-2 Project. It is obtained from inversion of Interferometric Synthetic Aperture Radar (InSAR) images acquired by the Sentinel-1 satellite. Although not reported herein for sake of brevity, sensitivity tests were carried out to check the impact of the Table 2 Crustal model used in numerical simulations. Adapted from Causse et al. (2021). ρ is the soil density, V P and V S are the P-and S-wave propagation velocities, respectively  Cornou et al. (2021) solution was eventually selected as reference source model because it was the one providing a slightly better match with recordings at the available stations. The used kinematic source parameters are summarized in Table 3, while the co-seismic slip distribution and the global Moment Rate Function (MRF) on the fault plane, in both time and frequency domain, are illustrated in Fig. 4. Locally, on each sub-fault, the MRF is defined according to Crempien and Archuleta (2015), with a rise time τ (defined as the local slip duration) equal to 0.5 s. The Fourier spectrum of the MRF turns out to be consistent with the Brune's model and shows a general consistency with the source spectrum adopted in Causse et al. (2021) (see Figure S3d). A sensitivity study has been performed to   (2021), a constant rupture velocity, V R = 1800 m/s, was adopted, corresponding to 85% of the V S of the top layer of the crustal model.
From a computational point of view, it is worth highlighting that the source modelling in SPEED takes advantage of a novel strategy, referred to as "not-honoring fault" (see Sangaraju et al. 2021, for the simulation of the 2016 Kumamoto earthquakes), specifically developed to account for finite-fault rupture models with arbitrarily complex geometries in a numerically efficient way. According to this approach, the mesh design does not need to incorporate the geometry of the fault plane (making the meshing operations timeconsuming and source-specific), but spectral nodes approaching the target fault rupture area are searched and loaded in order to reproduce the total seismic moment of the event to be simulated. Figure 5 illustrates the 3D spectral element numerical model, with indication of the finite-fault source area and the surface marking the boundary between the basin sediments and the underlying bedrock. The numerical domain extends over a volume of 45 km × 70 km × 8.5 km and it is discretized using a structured conforming hexahedral mesh with average length of the spectral elements ranging from about 120 m, at ground surface, to 550 m, at the bottom of the model. The numerical model of the valley has been generated through the not-honoring approach, where material properties across different geological sub-domains are given node by node at run time (Pelties et al. 2010). Referring to Sect. 4.2 for quantitative tests on the accuracy of the numerical solutions in the high-frequency range, the mesh was found to propagate accurately frequencies up to about 8 Hz. Using a fourth-order spectral polynomial degree (SD = 4), the total number of spectral nodes amounts to more than 80 Millions. Due to the large number of degrees of freedom, numerical simulations were performed on the Marconi 100 Cluster at CINECA, the largest high-performance computing center in Italy (www. cineca. it). The walltime for each numerical simulation is around 3 h on 128 cores of the Marconi 100 cluster, leading to a computational cost of almost 400 cores-hours.  Table 4 provides an overview of the numerical simulations that were performed in this study with the aim of (i) verifying and validating the numerical simulation and (ii) testing the impact of the 3D soil model on ground motion. With reference to (i), a simpler numerical mesh, with flat topography and crustal 1D layered structure was built.

Verification analyses with Hisada code
As a preliminary step of the modelling procedure, simulations were performed using the Hisada code, based on the analytical integration of Green's functions (Hisada and Bielak 2003). This code allows to compute the ground motions in a horizontally layered halfspace originating from a finite-fault kinematic source model.
Hisada's solution has been used in a preliminary phase to calibrate and validate the crustal model profile, the assumptions on the quality factor, the slip model on the . The simplified numerical model, without the basin shape and with the 1D crustal layering and flat topography, has been used for these tests (1D-pt model of Table 4). Figure 6 shows the recorded and simulated velocity time histories and corresponding Fourier amplitude spectra (FAS) obtained from SPEED and from Hisada's approach at station CRU1 and at a virtual receiver located at about 1 km from the source. All time histories have been low-passed filtered at 3 Hz, consistently with the frequency limit of Hisada solution. Herein a simple point-source model is adopted with an exponential source time function (rise time τ = 1.2 s), for consistency with the built-in source functions available in Hisada code. An excellent agreement between the two simulation techniques is found, especially in the near source, proving the accuracy of the numerical mesh. Agreement with recorded time series is good as well, both in amplitude, frequency content and arrival times. A detailed discussion on the misfit between simulations and recordings, in a broader frequency range, will be addressed in Sect. 5.

Numerical accuracy in the high-frequency range
Among the various approaches for numerical integration of the linear-elastodynamic wave equations, the spectral element approach enjoys a high accuracy, referred to as spectral accuracy, that was estimated to ensure an accurate wave propagation with slightly more than the Nyquist limit of 2 points per minimum wavelength (ppmw) for homogeneous soil conditions, up to about 4 ppmw in strongly heterogeneous materials (Faccioli et al. 1997). These estimates were based on verification tests on closed-form and/or reference solutions from literature. For a practical application, a proper check of the number of ppmw should be made for the specific case study, depending on the desired accuracy. For this purpose, Fig. 6 Verification of SPEED simulation (1D-pt model in Table 4 i.e, point-source 1D model, with flat topography and horizontal layers) against Hisada results. Simulated (blue: Hisada; red: SPEED) and recorded (where available) velocity time histories and corresponding Fourier Amplitude spectra for EW component at CRU1 station (R epi = 15 km, top) and a near field receiver (R epi = 1 km, bottom). All time histories are low-pass filtered at 3 Hz a convergence test was performed considering the numerical model described in Sect. 3.3 (3D-C21 model of Table 4), where, with the same discretization in terms of spectral elements, the spectral degree (SD) of each element was increased from SD = 1 (i.e., no internal Legendre-Gauss-Lobatto (LGL) node is present along each edge of the spectral element) up to SD = 5 (i.e., six LGL nodes within each side of the spectral element). In this way, the accuracy of the solution for SDj (j = 1, 2…5) can be checked by verifying at which frequency it departs significantly from the solution obtained with SDj + 1. Results of this test are illustrated in Fig. 7, where the acceleration waveforms and corresponding FAS are shown for selected stations on outcropping bedrock and in the basin for varying SDs, for the EW component of motion. This comparison indicates that, taking as a reference SD5, the solution with SD4 keeps close to SD5 up to about 7.5 Hz on outcropping bedrock and up to about 5 Hz in the basin. These should be considered as the reference accuracy limits However, it should also be pointed out that neither the input slip function nor the numerical model are detailed enough at high frequencies, which are dominated by smallscale effects mainly of stochastic nature. As it will be shown by comparing numerical results with records (see Sect. 5.2), the high-frequency decaying trend of simulated Fourier spectra is consistent with that of records. Because of such good agreement, we discarded the option to produce broadband results by coupling the low-frequencies from PBS to the high-frequencies produced by either stochastic methods or by Artificial Neural Network (ANN), such as proposed by Paolucci et al. (2018). Indeed, such hybrid approaches may not be theoretically well constrained for very shallow events, as it is the case of Le Teil earthquake. Moreover, in the case of ANN, a sufficient amount of records is necessary for training, which is not available for such shallow focal configurations. For this reason, we considered more physically sound to rely on the numerical content of the signal up to 8 Hz (i.e., signals were LP filtered below 10 Hz). Indeed, as shown in Fig. 7, both rock and basin receivers exhibit limited dispersion effects: the acceleration waveforms obtained for SD4 and SD5 present negligible differences in terms of peak values, phases and duration, with effects visible mainly in the coda of the signals (see zoom in bottom right panel of Fig. 7).

Overview of simulated results and comparison with records
In this section a summary of the reference simulation is given, comparing results from the 3D-C21 model of Table 4, with recordings and empirical GMMs, over the whole numerical domain.

Velocity motion and ground shaking maps
In Fig. 8, snapshots of horizontal ground velocity in the EW direction are shown, illustrating the patterns of seismic wave propagation and its interaction with the Rhône Valley. Basin induced amplification is noticeable both near the source, in the shallower portion of the basin, to the East of Montélimar (see snapshot at 7 s), and in the deeper Southern portion of the basin, such as near the OGLP station (see snapshots at 11 s). Figure 9 shows shaking maps of peak ground displacement (PGD), velocity (PGV) and acceleration (PGA) for vertical and horizontal components rotated along the strike fault normal (FN) and fault parallel (FP) directions. Recordings are shown as well, at reference stations, using the same palette. As discussed previously, to preserve the numerical results up to 8 Hz and allow consistent comparisons, simulated and recorded ground motions are filtered with a low-pass filter at 10 Hz. The maps of Fig. 9 point out an intense ground shaking in the immediate vicinity of the main rupture area, both in horizontal and vertical direction. The maps provide also a clear picture of the radiation pattern associated with a reverse focal mechanism and its interaction with complex subsurface geology. Updip directivity effects are visible on the hanging wall of the fault, yielding to a significant increase of ground motion amplitude, with maximum PGV and PGD values up to 0.4 m/s and 0.1 m, respectively. Above the fault projection, the PGAs simulated in this work range from 0.26 to 1.7 m/s 2 in the horizontal direction and from 0.45 to 2.6 m/s 2 in the vertical direction. These values are lower than the maximum PGAs found by Causse et al. (2021), up to 7.5 m/s 2 and larger than 1 g for the horizontal and vertical components respectively, probably because of the insufficient range of frequencies encompassed by our simulations and of the lack of high-frequency details of the slip distribution model.
Two prevailing directions of polarization of maximum amplitudes are noted, at azimuths of around 45° and 270° measured clockwise from North, because of the influence of the two shallow slip asperities (see Fig. 4) on rupture propagation across the fault. The influence of the basin sediments is clearly noticeable, especially in the PGA map, although limited in amplitude, because of the limited sediment depths in most portions of the basin except at sites around OGLP station. Basin effects are found also on the vertical motion, most likely because of the generation of surface waves at the basin edges. Vertical motions are high especially above the main slip area, with peak amplitudes which are comparable or even larger than the horizontal ones, consistently with observational evidences of large vertical-to-horizontal ratios in near-source conditions (Ramadan et al. 2021). FN and FP components show comparable amplitudes, although FN components tend to be larger than the FP ones on the surface projection of the fault rupture area.
From a qualitative point of view, the maps indicate a realistic spatial correlation structure of peak ground motion values, as it reflects the main features of the seismic wave propagation process (rupture and basin amplification effects). As expected, highfrequency measures of ground motion, like PGA, shows a spatial pattern with a smaller scale variability than PGV and PGD, the latter being associated with intermediate and long frequency components, respectively. Geostatistical analysis of the peak ground motion values through the semivariogram computation (not addressed herein) would confirm that the ground motion spatial correlation is consistent with observational evidences (for further details, see Infantino et al. 2021). Agreement with peak values of recorded motion is considerable, mostly in the horizontal directions.
Although not shown herein for sake of brevity, the simulated permanent vertical displacement on the surface projection of the fault plane (see, as a proxy for permanent ground uplift, the PGD-UD map in Fig. 9, bottom-right) reaches maximum values of  Figure 10 and Fig. 11 compare recorded ground motion with SPEED simulations, at the four stations: CRU1, on the basin edge to the North-East of the epicenter, and OGLP, inside the basin (Fig. 10), and ADHE and TRI2, on outcropping bedrock to the South-East of the epicenter (Fig. 11). Comparisons are shown in the time and spectral domain. For all stations, simulations are in satisfactory agreement with recordings in terms of amplitudes, arrival times and duration of motion. The vertical component tends to be in general more amplified, especially on later arrivals, both on bedrock and soil. At outcropping bedrock stations, ADHE and TRI2, simulations show on average a lower agreement with respect to CRU1 station; being on the other side of the basin, at larger distances from the source, most likely they are more influenced by uncertainties related to the adopted velocity model at regional scale (basin and crustal layers) which may be oversimplified. With reference to the crustal model, it is worth noticing that recent studies recommend the adoption of a laterally heterogeneous crustal model for the study area, with different properties in the eastern portion of the region (see e.g. Delouis et al. 2021;Bollinger et al. 2021). Inside the basin (at OGLP station), simulations are in very good agreement with observations for the main phases of ground motion, even though, in the horizontal direction, recordings exhibit reverberations with stronger amplitudes, not captured by the simulation, due to the simplified assumptions on the basin model. It is interesting to note that the simulations turn out to reproduce relatively well the high energy content of vertical motion at periods larger than about 0.5 s (see e.g. CRU1 spectra). Although simulations tend to overestimate recorded vertical motion at long periods, they can capture the local peak in the response spectra  (Ramadan et al. 2022) confirm that, for Le Teil earthquake, vertical-to-horizontal ratios of response spectral ordinates tend to be systematically larger than those predicted by empirical GMMs in the range of intermediate to long periods.

Goodness-of-Fit Scores
The overall performance of the numerical simulations was quantitatively estimated through the Goodness-of-Fit (GoF) criteria proposed by Anderson (2004), considering ground motion parameters of interest for earthquake risk applications, namely: PGA, PGV, PGD, Acceleration Response Spectra (SA) at selected periods (T = 0.2, 0.5, 1 and 2 s) and Cumulative Absolute Velocity (CAV). The scores were computed in the frequency range up to 8 Hz. Figure 12 (left, model 3D-C21) shows the individual scores associated with the aforementioned parameters for the set of four reference stations, for both horizontal geometric  Fig. 10 a, for ADHE (R epi = 18 km) and TRI2 (R epi = 24 km) stations 1 3 mean (GMH) and vertical (UD) components. An overall good-to-excellent agreement is found for the horizontal components for all stations and all parameters, while for the UD component, slightly worse scores are found, due to the higher uncertainties in constraining the P wave propagation model. The closest station (CRU1) shows the best performance on the vertical component.

Sensitivity analyses on rise time
In order to verify the accuracy of the results in the high frequency range, especially in the vicinity of the fault, we have repeated the numerical simulations by varying the rise time of local slip time function from the reference value τ = 0.5 s (3D-C21 model of Table 4) to 0.3 s (3D-C21-0.3 model). Figure 13 shows the results of 3D-C21 and 3D-C21-0.3 models (in red and blue, respectively) at stations CRU1, OGLP and TRI2, in comparison with records (in black), in time and frequency domain, for EW component of motion. It is clear that the decrease of the rise time leads to an increase of high-frequency energy radiation especially at frequencies larger than 1 Hz. However, such an increase of high-frequency  Table 4) EW velocity time histories and corresponding Fourier spectra at CRU1, OGLP and TRI stations. All time histories are low-pass filtered at 10 Hz harmonics is in excess with respect to observations, with an overall decay of the performance of simulation. The lower performance of the 3D-C21-0.3 model is demonstrated by the GoF analysis shown in Fig. 12 of previous section, where GoF scores from 3D-C21 (left) and 3D-C21-0.3 (right) models are compared. It is clear that, for lower rise time, the GoF scores, evaluated in the frequency range up to 8 Hz, show a strong reduction for all stations and for all metrics under consideration. For this reason, model 3D-C21 is selected as reference model in our study, although, as mentioned previously, near-fault PGA values tend to be smaller than the ones found by Causse et al. (2021) from forward numerical simulations of different realizations of slip distribution along the fault, supported by in-situ observations of funeral slab displacements. Figure 14 shows horizontal PGA, SA(0.2 s), SA(0.5 s) and SA(1 s) as a function of Joyner-Boore distance, R jb , from 3D-C21 SPEED simulations (black: on outcropping bedrock; grey: inside the basin) as well as from recordings (green dots) and from the GMM by Kotha et al. (2020), K20, for shallow crustal events (violet: rock, with V S30 = 2100 m/s; red: basin, with V S30 = 500 m/s), in its ergodic formulation. V S30 values of empirical predictions are selected to be consistent with the Vs profiles implemented in the numerical model. For both simulations and recordings, the median (RotD50) values of spectral accelerations over all orientations (Boore 2010) are considered for consistency. The same comparison is shown in Fig. 15 but for the vertical component, using the GMM by Stewart -Numerical simulations are suitable to fill in the gap left by the recordings in the proximity of the seismic source: while the records cover, with a very limited sampling, only distances beyond about 15 km, numerical simulations provide a detailed picture of ground motion in the near-source region (0-15 km), with relevant implications for constraining the seismic input for the NPPs; -In general, in the distance range between 15 and 25 km, a satisfactory agreement is found between simulated and recorded spectral values, although some discrepancies are found at CRU1 station, especially for the horizontal PGA and vertical SA(0.5 s); -Simulations suggest a tendency of near-source ground motions to be lower than the GMMs model at short periods (particularly for PGA), while a reverse trend is found at long periods (≥ 1.0 s), where simulated ground motions exceed the empirical predictions. Such a trend is consistent with the findings from residual analysis of Le Teil recordings with respect to recorded datasets in France (Ramadan et al. 2022). Furthermore, it should be noted that PGA may be underestimated because of the frequency limit (8 Hz) of the PBS; -The systematic underestimation of the empirical GMMs for the long-period spectral ordinates is even more pronounced for the vertical components of ground motion, for which both simulations and recordings tend to be above the 84 th percentile of the empirical predictions. As mentioned previously, this feature was found also in independent empirical studies on the French recorded dataset, indicating that, for the Le Teil earthquake, long periods (say at around 1 s) vertical spectral accelerations

Basin amplification effects
One of the main advantages of PBS is that, unlike for records, site conditions are completely known. While the complexity of the small-scale variability of ground properties cannot be portrayed in detail, unless spatially correlated stochastic fields are applied to the average values of wave propagation velocities (see Paolucci et al. 2021, for an application to the PBS of ground motions from induced seismicity in Groningen, Netherlands), PBS are suitable to investigate site amplification effects from a variety of viewpoints, such as: (i) the variability of site amplification with respect to different outcropping bedrock stations; (ii) the comparison with the results from 1 and 2D simulations (see e.g. Smerzini et al. 2011, for an application to the Gubbio basin, Italy); (iii) their repeatability and scenario dependence both in the linear and non-linear ranges (see e.g., Stafford et al. 2017).  Table 4 Although SPEED allows for consideration of a relatively simple, albeit effective, non-linear visco-elastic model (Stupazzini et al. 2009), due to the relatively low levels of seismic excitation we will consider in this section only linear visco-elastic site amplification effects. To this end, the focus is on the cross-section shown in Fig. 16 (the azimuth of the cross-section is the same as the fault strike), where velocity time histories at selected locations are illustrated, on the top panel, both in the basin and at outcropping bedrock, and, on the bottom, at all receivers along the cross-section (the EW component was chosen as reference for these investigations, for simplicity). The latter plot allows to highlight the complexity, as well as the 3D nature, of seismic wave propagation in the basin. As clarified also from Fig. 8, the western portion of the cross-section is the one that first experiences ground motion due to the lower source-to-site distance, but afterwards the presence of the basin increases the complexity of the overall seismic wave field, with prominent amplification effects towards the basin center.
To quantify the basin-induced amplification, regardless of the reference station, a second PBS, with the same reference kinematic source model (see Sect. 3.2), was carried out but without the presence of the basin. In such simulation, referred to as "3D-C21-R" (where "R" means "rock") in Table 4, the dynamic properties of the basin are replaced by those of the outcropping bedrock. In Fig. 17 (left), the acceleration time histories at D and C receivers (shown in Fig. 16) are plotted, clearly pointing out the increase in amplitude and in duration of ground motion as well as the elongated dominant period with respect to the case without basin. Such prominent amplification is clearly shown in terms of the corresponding acceleration response spectra, highlighting the significant long period amplification especially at the center of the basin.
The quantification of site effects is further explored in Fig. 18, where the Fourier spectral ratios (FSR) and Response spectral ratios (RSR) are considered, with reference to the receivers C and D: (a) label "3D" refers to the spectral ratios obtained by dividing the results of the simulation 3D-C21 over the 3D-C21-R; (b) label "1D" refers to the 1D theoretical amplification function with the local stratigraphy below the corresponding receiver; for the RSR, the accelerograms at C and D were computed by 1D convolution using as Fig. 17 Comparison of EW acceleration time histories with and without basin (3D-C21 and 3D-C21-R simulations, respectively) and corresponding acceleration response spectra (SA), for the D (top, blue lines) and C (bottom, red lines) receivers of Fig. 16. Basin thickness ranges from 250 m (receiver C) to 700 m (receiver D) while V S30 = 500 m/s at both sites input motions the corresponding accelerograms computed in the without basin case (3D-C21-R); (c) labels "C/B" and "C/F" (and similarly for receiver D) refer to the spectral ratios computed from 3D-C21 run with respect to reference stations B and F, located on the left and right side of the basin respectively (see Fig. 16).
Several remarks can be made based on the inspection of Fig. 18: -with some exceptions in relatively small frequency intervals, the FSR show that there is an overall good agreement between the 1D and 3D amplification functions, indicating a moderate influence of 3D site effects; -in agreement with studies aiming at quantifying the aggravation factors on the response spectra related to complex 2D/3D geological configurations (e.g., Chávez-García and Faccioli 2000; Riga et al. 2016), the 1D solution tends to overestimate the 3D one close to the basin edge (receiver C), while the opposite occurs at the basin centre (receiver D). However, comparison of RSR at short periods suggests that PGA from 3D simulations tends to be smaller than in the 1D case; -if site amplification functions are computed with respect to a reference station, the location of such station with respect to the basin is critical (as also discussed in Smerzini et al. 2011): namely, spectral ratios with respect to receiver F show some sharp anomalies, since F lies on the other side of the Rhône basin with respect to the source, so that a part of the frequency content, including low and high frequencies, is filtered out by Fig. 18 Comparison of amplification functions computed from Fourier spectral ratios (top) and Response spectral ratios (bottom) at the selected receivers C (left) and D (right). The plots show the spectral ratios from 3D-C21 simulations considering B and F as the reference stations on outcropping bedrock (see Fig. 14 for their location). 1D refers to the local 1D theoretical amplification functions, while 3D are the spectral ratios obtained by dividing the results of the 3D-C21 run (with basin) to the 3D-C21-R one (without basin). For the response spectral ratios, the accelerograms at C and D were computed by 1D convolution using as input motions the corresponding accelerograms computed in the without basin simulation the presence of the basin; spectral ratios with respect to receiver B are instead closer to the 1D and 3D solutions.
More generally these results confirm that, when complex geological configurations lie in the vicinity of active faults, the main features of seismic response cannot be reliably captured by standard approaches owing to the variability of the source-to-site ray paths affecting wave propagation. In these conditions, especially in case of critical structures such as NPPs, PBS seem to be an effective way to predict the regional as well as the site-specific features of the seismic response. It is advisable that temporary earthquake networks are installed in the region to provide further observational evidences on local site amplification in the Rhone valley, to be used also for verification and validation of the 3D numerical model.
Further investigations are planned, starting from this case study, aiming at evaluating the variability of site effects from different realization of earthquakes from the same source, with variable magnitude and slip distribution, and from other sources in the investigated area, with different distance and azimuth from the site.

Conclusions
This paper presents a case study of validated physics-based numerical simulation (PBS) at regional scale. The target area is an industrial region in South-Eastern France within the Rhône Valley, which hosts several operating nuclear power plants (Cruas and Tricastin). This low-to-moderate seismicity area was hit by an unusually shallow (1 km depth) M W 4.9 earthquake rupturing the La Rouvière fault, near Le Teil, at about 15 km distance from the Cruas NPP. The recordings available for the Le Teil earthquake, although limited in number, were used for validation of the PBS. The numerical code SPEED was used for this purpose. In the recent past, this numerical code underwent several successful validations with earthquakes in a wide range of magnitudes, spanning from 3 to 6.5, and in different countries worldwide.
A model of the Earth's crust was constructed, including the La Rouvière fault and the lower Rhône Valley, with a size of 70 km × 45 km × 8.5 km and a total of about 82 millions nodes using a spectral degree SD = 4. A convergence test on the accuracy of the numerical solution for varying SDs indicated the presence of moderate dispersion effects in a broad frequency range up to around 8 Hz, which was regarded as the maximum reliable frequency of PBS. For this reason, the numerical solutions were eventually filtered below 10 Hz and no hybrid technique was applied to produce broadband signals.
When comparing simulations with records by means of quantitative GoF metrics, a good-to-excellent agreement was found in a frequency range up to 8 Hz for all stations and for several horizontal ground motion parameters. Although some limitations are still present in simulations owing to simplifying assumptions regarding the relatively lowfrequency source model as well as the basin model, the set of sanity checks performed in this study, ranging from the GoF analysis, attenuation with distance to ground shaking maps and site amplification features, indicates that 3D PBS may provide realistic broadband predictions of earthquake ground motion at regional scale. On the other hand, a thorough analysis on ground motion prediction in the immediate vicinity of the fault, as addressed by Causse et al. (2021), is beyond the scope of this work.
This work also demonstrates that, with limited recordings available (as for the low-tomoderate seismicity region of Montélimar) and even without very detailed 3D velocity and source models, PBS can be employed to shed light on a variety of aspects related to ground motion modelling, poorly addressed by ergodic empirical GMMs, spanning from the prediction of ground shaking intensity and spatial variability in the proximity of the seismic source, to region-and site-specific features of ground response. The investigation of such aspects is particularly relevant when performing seismic risk evaluation for strategic and critical structures, infrastructures and industrial plants, such as NPPs, the failure of which during an earthquake may endanger safety of population and cause environmental disasters. For such facilities, the definition of seismic hazard at long return periods may take advantage of the capability of 3D PBS of accounting for region-and site-specific features of ground shaking, especially in near-source conditions, that are typically not accounted for by ergodic empirical models.
In addition to these advantages, validated PBS act as a numerical laboratory, where realistic time histories of ground motion can be obtained, under fully controlled knowledge of the dynamic parameters affecting the seismic wave propagation and of the slip distribution along the fault. Such results may further provide a wealth of information on key issues related to earthquake ground motion, such as explored in this paper with the investigation on the basin-induced site amplification within the Rhȏne Valley, that was shown to provide physically sound results in a broad frequency range. Taking advantage of on-going studies on the collection and interpretation of geological and geophysical data of the Rhȏne Valley, the 3D PBS may be improved in future works and the prediction capability of a refined velocity model, with respect to a simplified one, may be quantified and discussed.