Optimization of a Simulation Code Coupling Extended Source (k−2) and Empirical Green’s Functions: Application to the Case of the Middle Durance Fault

We developed a ground-motion simulation code base on extended rupture modeling combined with the use of empirical Green’s functions (EGFs), adapted for low-to-moderate seismicity regions (with a limited set of EGFs), and extended its range of applicability to the lowest source-to-site distances. This code is based on a kinematic source description of an extended fault and is designed to allow complex fault geometries and to generate a ground motion variability in agreement with that of the recorded databases. The code is developed to work with a sparse set of EGFs. Each available EGF is therefore used in several positions on the rupture area. To be used in positions different of their original position, we applied to the EGFs some adjustments. In addition to the classical adjustments (i.e. time delay correction, geometrical spreading correction and anelastic attenuation correction), we propose here a radiation pattern adjustment. The effectiveness of it is tested in a numerical application. We showed noticeable improvements at the lowest distances, and some limitations when approaching the nodal planes of the subevents the recording of which were used as EGFs. We took advantage of the development of this code, its ability to work with a sparse set of EGFs, its ability to take into account complex fault geometries and its ability to master the general variability, to perform a ground-motion simulation scenario on the Middle Durance Fault (MDF). We perform simulations for a hard rock site (VS30 = 1800 m/s) and a sediment site (VS30 = 440 m/s) of the CEA Nuclear Research Site of Cadarache (France), and compared the computed ground motion with several ground motion prediction equations (GMPEs). The GMPEs slightly underestimate the sediment site but strongly overestimate the ground motion amplitude on the hard rock site, even when using a specific correction factor which adapts GMPEs predictions from rock site to hard rock site. This general ascertainment confirms the need to continue efforts towards the establishment of consistent GMPEs applicable to hard-rock conditions.


Introduction
The aim of this study was to develop groundmotion simulation code based on extended rupture modeling adapted for low-to-moderate seismicity regions, and to extend its range of applicability to the lowest source-to-site distances (i.e. typically distances smaller than the geometric extension of the fault, for which the variation of azimuth between the site and the different extremities of the rupture are non-negligible). This code is based on a kinematic source description of an extended fault according to the k -2 model (Causse et al. 2009;Dujardin et al. 2016;Del Gaudio et al. 2018), which is combined with the use of empirical Green's functions (EGFs). It is freely available on github: https://github.com/ adujard/K2_FMD/.
Following the EGF approach (Hartzell 1978), small magnitude events are used as EGFs to allow implicit consideration of the characteristics of the propagation medium (e.g., anelastic attenuation, geometrical spreading, site properties). The advantage of this approach is the possibility to simulate ground motion in a wide frequency range when purely numerical approaches are limited by partial Electronic supplementary material The online version of this article (https://doi.org/10.1007/s00024-019-02309-x) contains supplementary material, which is available to authorized users.
knowledge of the small-scale heterogeneities of the propagation medium. However, the deployment of this method in low-to-moderate seismicity areas (as for metropolitan France) is often limited by the number of available small-magnitude events that are usable as EGFs (limited, in addition to the quality of the recording, by range of magnitude, position, depth, focal mechanism; see Dujardin et al. 2016, for details).
The simulation code is then essentially developed for applications with a limited and sparse set of EGFs. This involves being able to use an EGF at a location different from its original position. Therefore, the classical adjustments have been implemented to correct for theoretical differences in geometrical spreading, anelastic attenuation, and travel-time differences between the initial position of the EGF and the position at which it is moved to, to satisfy the lack of data. In addition to these classical adjustments, we implemented a radiation pattern correction, which constitute the most important adjustment and the major contribution of this study.
This adjustment is all the more important because the radiation pattern has a major influence on groundmotion variability at short distances (Dujardin et al. 2018). This correction also makes it possible to improve the prediction at the lowest source-to-site distances, where the applicability of ground-motion prediction equations (GMPEs) is limited due to the lack of data at short distances. Indeed, several studies have highlighted that the use of GMPEs outside the magnitude and distance ranges covered by the datasets from which they are derived leads to missestimation of the ground motion (Bommer et al. 2007;Chiou et al. 2010;Douglas and Jousset 2011). Moreover, close to the fault, the ground motion amplitude is affected by the saturation effect (e.g., Yenier and Atkinson 2014). Although its influence is taken into account through most of the GMPE functional forms (see Fukushima et al. 2003;Zhao et al. 2006), it is only an average consideration for a set of different events and different azimuths. This approach is therefore not suitable for specific site prediction, and does not allow the large azimuth variability of the ground motion and the saturation effect to be taken into account (see Dujardin et al. 2018). Furthermore, the radiation pattern correction allows some complex fault geometries to be taken into account, with variations of strike and dip angles. This is the case for the geometry of the Middle Durance Fault (MDF) (Guyonnet-Benaize et al. 2015), which constitutes our application case here (see Sect. 4).
Finally, this simulation code is developed in the framework of blind ground-motion prediction in lowto-moderate seismicity areas. This means that it is designed for simulation of a hypothetical future event for which the source characteristics (e.g., stress drop, rupture velocity, rupture dimensions) are a priori unknown. Therefore, the input rupture parameters are defined by the probability density functions, the variability of which can be constrained to reproduce the ground-motion variability observed (for the target event characteristics) in strong-motion databases.
The last part of this paper presents an application of the simulation code, to generate a suite of groundmotion time series for potential ruptures on the MDF, which is located in southeastern metropolitan France, using stations located at the CEA Nuclear Research Site of Cadarache (France). The MDF is a complex system of several segments that has a global extension of about 60 km (Volant et al. 2000), and is mainly oriented North-East to South-West. This application was chosen for several reasons. On the one hand, the fault passes close to the CEA (* 5 km), and on the other hand, this fault has real seismogenic potential that has been expressed in the recent past by several earthquakes (Fig. 1) of magnitude around 5 (Lambert et al. 1998). Predicting reliable ground motion at or close to the Cadarache site associated to the occurrence of hypothetical earthquakes on the MDF constitutes a scientific challenge, and is clearly an issue for the design and assessment of future or existing facilities in the context of nuclear safety (Berge-Thierry et al. 2017a).
case, there are no available recordings and the physical characteristics of the source are unknown, except for the seismic moment M 0 and the focal mechanism, which are postulated. We first define the size of the fault on which the rupture is expected. Then, the size of the rupture area on this fault, which is assumed to be rectangular, is automatically calculated from the stress drop (Dr) and the seismic moment M 0 , as originally proposed by Herrero and Bernard (1994).
The stress drop needs to be chosen by the user, for several reasons. First, although the stress drop has been considered as constant since the fundamental work of Aki (1967), it might vary from one region to another (Chouet et al. 1978). Secondly, the variability of the stress drop is significant, with variations generally between 0.1 and 100 MPa (Allmann and Shearer 2009). Moreover, the hypothesis of earthquake self-similarity is questionable, and it is suggested that the constant stress-drop condition is not appropriate for a large range of magnitudes (Mayeda and Walter 1996;Beeler et al. 2003;Kanamori and Rivera 2004;Drouet et al. 2011).
The dimensions of the rupture area are derived as follow: from the input stress drop (Dr), the theoretical corner frequency (fc) is derived following Brune (1970): where V S is the S-wave velocity (m/s), M 0 is the seismic moment (Nm), the stress drop is expressed in Pa, and k = 0.37 (Brune 1970). Then, according to the following approximation (Hanks 1979;Hanks and McGuire 1981): is the rupture duration and V R is the rupture velocity (m/s), we derived the length (L) and width (W) of the rupture area by assuming that the characteristic size of the rupture is . Thus, only the ratio between L and W is necessary to derive the dimensions of the rupture area. V R depends on V S in the vicinity of the fault, and it commonly varies between 0.7* V S and 0.85* V S (Heaton 1990). V S in the vicinity of the fault is also used to derive the differences in travel-times between the different parts of the rupture area and the target station. Both V S and the ratio between V S and V R are parameters to be chosen by the user.

Static Slip Generation
Once the rupture dimensions are defined, the static slip distributions of the source are generated in two steps, as the low and high frequency parts of the static slip are constrained separately. The low frequency part is set to a constant value over the rupture area, as the mean slip, which is derived from the relation between the seismic moment and the rupture dimensions: where A ¼ L Â W is the rupture area, l ¼ qV 2 S (q = 2.7 g/cm 3 ), and V S is the S-wave velocity. The high frequency part of the static slip distributions (representing slip heterogeneities) is defined in the wavenumber domain for any wavenum- . Following Herrero and Bernard (1994), this high frequency part should have an asymptotic decay in k -2 beyond the corner wavenumber k c . The high frequency part of the displacement spectrum is defined for the rectangular , the CEA Cadarache center, and the main historic seismic events (Lambert et al. 1998) Vol. 177, (2020 Optimization of a Simulation Code Coupling Extended Source (k -2 ) 2257 fault plane similarly to Somerville et al. (1999) and Gallovič and Brokešová (2003): where k x and k y are the wavenumbers in the strike and the dip directions, respectively, U k x ; k y À Á is the phase spectrum, which is randomly defined, and kc x and kc y are the corner wavenumbers in the strike and dip directions. In our code, kc x and kc y can be used and interpreted in two ways. First, these corner wavenumbers can be chosen by the user to constrain the characteristic dimensions of slip heterogeneities along the strike and dip (D cara x;y ¼ 1 kc x;y ). This approach therefore makes it possible to control the roughness degree of the static slip (see Causse et al. 2009).
In the second interpretation, the corner wavenumbers associated to the strike and the dip directions are defined as kc is the characteristic rupture dimension, and T RUP is the rupture duration, which is proportional to the corner frequency: f C ¼ 1 T RUP (Hanks 1979;Hanks and McGuire 1981). This leads to the following relationship between k C , the theoretical corner frequency (f C ), and the rupture velocity (V R ): Following this second approach, the corner frequency of the simulated event (with one given moment magnitude, M 0 ) does not depend on the rupture dimensions (i.e., L, W) nor on the rupture velocity (V R ), but only on the static stress drop chosen as an input parameter. Thus, in this approach, the corner frequency variability is directly controlled by the variability of the input parameter Dr.

Spatial Sampling
The rupture area is discretized into subfaults where their sizes are defined according to the target maximum frequency of the simulations. Indeed, the static displacement amplitude spectrum (Eq. 3) is valid only up to the Nyquist wavenumber (k NY ), which depends on the sampling rate in the spatial domain (subfault size SF dim ): The generated ground-motion spectra are then limited to the maximum frequency defined by f kmax % k NY Â V R . The subfault size SF dim is then automatically defined in the code according to the maximum target frequency (f kmax ): Note that our approach differs from the EGF formulation based on scaling laws between large and small earthquakes, in which the subfault size depends on the EGF seismic moment (e.g., Irikura and Kamae 1994;Miyake et al. 2003;Asano 2018). As explained in Sect. 2.5, the EGFs are processed so that they can be treated as approximations of ''real'' Green's functions (flat displacement amplitude spectrum and unit seismic moment).

Rupture Kinematics
The rupture kinematics describes the space-time evolution of slip on the fault plane. First, the rupture time at position x; y ð Þ on the fault is defined as T x;y = D nuc x; y ð Þ=V R , where D nuc x; y ð Þ is the distance to the nucleation point and V R is the average rupture velocity. Following the approach of Hisada (2000Hisada ( , 2001, we introduce rupture velocity variations through perturbations of the rupture time, distributed over the rupture area, also according to a k -2 model. The rupture time is given by: where DT R x; y ð Þ is the rupture time variation, expressed as a percentage where its maximum controls the roughness level of the propagation front. Following the Hisada (2001) formulation, the amplitude spectrum of the rupture time variations on the rupture area is given by: where kc Tx and kc Ty are the corner wavenumbers in the strike and dip directions, respectively. As for the high frequency slip distribution (Eq. 3), the corner wavenumbers control the characteristic size (Dim x;y DT R ) of the time perturbation in the strike and dip directions: Second, the slip rate function is defined as the sum of the isosceles triangles as proposed by Hisada (2001). The slip rate function can be parameterized by three parameters : the slip rate function duration s rise (defined as s max in Hisada (2001)), the number of summed triangles (Nv) and Ar which corresponds to the ratio of the area of the j?1th triangle with respect to the ratio of the jth triangle (i.e. Ar = A j?1 /A j ). In the following of the paper we use Nv = 4 and Ar ¼ ffiffi ffi 2 p for a M 6.0 event. Hisada (2001) showed that is has two characteristics frequencies: f 1¼ 1= 2s rise ð Þ, where s rise is the rise time (i.e., the slip rate function duration) and f max¼ 1=s 1 , where s 1 is the duration of the first triangle. s rise is supposed to be constant over the rupture area, and is defined according to Somerville et al. (1999): where M 0 is in dyne.cm. Finally, we obtain the absolute source time function or moment rate function M 0 in the time and frequency domains by summing up the contributions of each subfault (each fault point): Or as described in the frequency domain (see Eq. 3 in Hisada (2001)): where x and y denote the along-strike and along-dip directions, respectively. l is the rigidity, D(x,y) is the static slip (see Sect. 2.2), T(x,y) is the rupture time, and SRF t ð Þ is the normalized slip rate function (that is corresponding to unite slip). Hisada (2001) shown empirically that such a slip rate function associated with the above-mentioned k -2 distributions of static slip and rupture time perturbations, results in a the classically observed x À2 decay of the moment rate function. More precisely, the slip velocity models, built by superposing equilateral triangles introduce characteristic frequencies: f 1¼ 1= 2s rise ð Þ and f max is imposed by the minimum duration among the triangles. Then, the amplitude spectra of the absolute source time function fall off x À1 before f 1 and as x À2 between f 1 and f max . Note that an additional characteristic frequency f kmax is introduced by the fault discretization. As explained in Sect. 2.3, the subfault size is chosen to reach the target frequency fixed by the user.

Time Series Generation and Green's Function Corrections
Following the discrete representation theorem (Aki and Richards 2002), the simulated acceleration U r * ; t for a station at position r * is expressed as: where R x; y; t ð Þ represents the contribution to the moment rate function at position (x,y), and FG x;y r; t ð Þ is the Green's function in acceleration associated to the same subfault. As mentioned above, there is no specific EGF associated to each subfault. The Green's function FG x;y r; t ð Þ at position (x,y) is then approximated using the closest available EGF, FG 0 r; t ð Þ. First, the EGF is deconvoluted from its source spectrum, so as to be considered a Green's function. This deconvolution is performed in the frequency domain, by preserving the EGF phase spectrum. The amplitude source spectrum S(m 0 , f c ) is approximated according to the x -2 model (Brune 1970), which depends only on the moment magnitude m 0 and the corner frequency f C of the EGF: Vol. 177, (2020) Optimization of a Simulation Code Coupling Extended Source (k -2 ) Secondly, we apply several adjustments to this Green's function when it is shifted from its original position. We first apply a time shift to the initial Green's function FG 0 r; t ð Þ, which is equal to the travel-time difference Dt(x,y) between its initial and new position to the considered station: with the travel-time difference Dt(x,y) defined as follows: where R SF x;y STA r ð Þ is the distance between the subfault at position (x,y) on the rupture area and the station considered, and R FG 0 STA r ð Þ is the distance between the initial Green's function and the station.
The difference in the geometrical spreading (defined as 1/R c ) is corrected as follows: Anelastic attenuation is also accounted for, with it modeled by the quality factor Q(f). The correction is applied in the frequency domain as follows: where FFT is the fast Fourier transform operation, iFFT is the inverse transformation, and Q f ð Þ ¼ Qsf a is the quality factor estimation for the S wave, where f is the frequency.
Finally, we implement a frequency dependent radiation pattern correction, as proposed by Pitarka et al. (2000). We use the Aki and Richards (2002) coefficients, computed for S-waves that propagate in a homogeneous infinite elastic medium in the farfield approximation (see Aki and Richards 2002, Eq. 4.33). The Green's function time series first needs to be expressed in the local coordinate system (r h U), specific for each station, which naturally brings out the radial and transverse components of motion (see Aki and Richards 2002, Fig. 4.4). First, the coordinates are rotated from the ENZ to the (y1 y2 y3) coordinate system which is centered on the fault (Fig. 2). The axes (y1 y2 y3) correspond to the axes (x2 x1 À x3) of the Aki and Richards (2002) coordinate system (see Aki and Richards 2002, Fig. 4.4). We consider that for a strike (U S ), a dip (d), and a rake (k) of 0, the y1 ¼ x2 ð Þ direction corresponds to the East, the y2 ¼ x1 ð Þdirection to the North, and the y3 ¼ À x3 ð Þdirection to the Z direction ( Fig. 2). Then the temporal series of Green's functions expressed in the (y1 y2 y3) coordinate system are given by: Representation of the (y1 y2 y3) coordinates on the ENZ components orientation. The (x1 x2 x3) coordinates used by Aki and Richards (2002) are also indicated. This case is represented for a strike, a dip, and a rake of zero. The position coordinate h and U used in Eq. (4.33) of Aki and Richards (2002) are also indicated, as are the orientation of the (r h U) local coordinate system.
Note that we also derive the coordinates of the station in the (y1 y2 y3) system coordinates by rotating the distance in the North, East and vertical directions between the source and the station: where STA E;N;Z and EGF E;N;Z are the positions of the station and the EGF in the ENZ system, respectively, and STA y1;y2;y3 is the position of the station in the (y1 y2 y3) coordinate system. In this way, it is possible to compute the spherical coordinates h and U of the station position and rotate the Green's function time series into the (r h U) local coordinate system ( Fig. 2; see also Aki and Richards 2002, Fig. 4.4): Once the Green's function is expressed in the (r h U) coordinate system, we can compute the farfield coefficients A FP ini and A FS ini of the P waves and S waves, respectively (see Aki and Richards 2002, Eq. 4.33), associated to the initial position of the EGF. As these are expressed in the (r h U) coordinate system (i.e., withr 1; 0; 0 ð Þ;ĥ 0; 1; 0 ð ÞandÛ 0; 0; 1 ð Þ), the A FP coefficient has only one non-zero component in ther direction. The A FS coefficient has two nonzero components: A SH ¼ A FS ÁÛ associated to the SH-waves, and A SV ¼ A FS Áĥ associated to the SVwaves. The coefficients A FP ij and A FS ij are then derived for the subfault position on which the Green's function is moved: U S;ij , d ij , and k ij for the strike, dip, and rake associated to the subfault i; j ð Þ, and h ij and U ij for the relative position in the (r h U) coordinate system of the station to the ij ð Þ subfault. Once expressed in the (r h U) coordinate system, the temporal series can be corrected. However, due to the isotropic nature of the high-frequency radiation pattern (Liu and Helmberger 1985;Takenaka et al. 2003;Takemura et al. 2009;Sawazaki et al. 2011;Kobayashi et al. 2015), a frequency-dependent correction is considered. This correction is thus defined in the frequency domain, as follows: where FFT and iFFT are the fast Fourier transform and the inverse operation, respectively, the dot is a term to term multiplication, and RPC r;h;U is the radiation pattern correction (Fig. 3) that is expressed as: where f is the frequency, and A r;h;U is the correction coefficient, which is defined as: Finally, a threshold for this correction is defined. If one of the EGF (i.e. FG 0 r; t ð Þ) radiation pattern amplitudes (A FP ini Ár, A SH ini , A SV ini ) is \ 10 -1 , no amplitude correction is performed. This limit is intended to avoid having to divide by a number closed to 0, and to avoid strong amplification of the noise when working with real data, which would pollute the simulations.
Then the temporal series are rotated back into the ENZ system coordinates, according to the strike, dip, and rake of the ij ð Þ subfault (U S;ij d ij and k ij ) and the position of the ij ð Þ subfault relative to the station (h ij and U ij ): A purely numerical test is now carried in Sect. 3. to test the relevance of this correction.

Description of the Experiment
To test the efficiency of the radiation pattern correction, the experiment of Dujardin et al. (2018) is reproduced here, in which the impact of the radiation pattern on the saturation effect of the ground-motion peak values is highlighted. According to Dujardin et al. (2018), a M 6.0 event is generated (i.e., M 0 = 1.1220.10 18 Nm). The fault is defined as a plane with a strike of 0°, and a dip of 60° (Fig. 4). The dimensions are fixed for every simulation to L = 12650 m and W = 7900 m, with fixed input parameters: Dr ¼ 1:0 MPa, f C ¼ 0:17 Hz, V S ¼ 3600 m/s, and V R ¼ 0:7 Â V S ¼ 2520 m/s. The subfault dimensions are fixed to SF dim ¼ 36 m , so that the ground-motion simulations are valid up to f kmax ¼ 35 Hz (Eq. 6).
A population of 20 static slip distributions is generated using the k -2 model (Eq. 3). The rupture times are disturbed according to the model proposed by Hisada (2001) (Eq. 8). The characteristic dimensions of the rupture time variations are fixed between 30% and 70% of the L and W dimensions of the rupture area (Eq. 9). The amplitudes of the rupture time variations are normalized, so as not to exceed 10%, to avoid a propagation front that is too rough, which would generate additional high frequencies beyond what is expected by the x -2 model. The nucleation is located at depth, at distances of 0.15*L along the strike direction, and 0.8*W along the dip direction (Fig. 5). To facilitate the reading of the results, no travel-time correction (Eqs. 15, 16) is considered for any of the cases presented here, to suppress any directivity effects (Dujardin et al. 2018). Simulations are performed for stations located to the North only, at azimuths from -90°to 90°with respect to the fault center, every 22.5° (Fig. 4). The stations are located at rupture distances R RUP of 3, 6, 10, 20, 40 and 70 km.
Three simulation cases are considered. In the first, one unique Green's function is considered, which is located at the center of the rupture area. When this one is shifted from its original position, only the geometrical spreading correction (Eq. 17) and the anelastic attenuation correction (Eq. 18) are considered. In the second case, a single Green's function is also used, but radiation pattern correction (Eqs. 19-25) is added. The third case represents the reference solution, and it allows the effectiveness of this radiation pattern correction to be judged. In this case, the radiation pattern complexity is reproduced by the use of 45 Green's functions distributed every 1.5 km along the strike and dip of the rupture area (Fig. 5).
The Green's functions are computed numerically using the discrete wavenumber method (Bouchon 1981). In every case, these Green's functions are computed for the focal mechanism defined by the strike, dip, and rake angles of 0°, 60°, and 0°, respectively. The propagation medium is assumed to be a homogeneous half space where V P = 5000 m/s and V S = 3600 m/s, and the anelastic attenuation defined through the quality factor is fixed at Q P = 50f 0.2 for P waves and Q S = 200f 0.3 for S waves. As the propagation medium is homogeneous, the radiation pattern correction is in this case applied to the whole frequency range (Eq. 23).
The results of the three cases are compared in terms of the peak ground acceleration (PGA) versus the R RUP distance. The PGA for each station is the mean PGA over the 20 rupture realizations. The time series are initially low-pass filtered, at under 35 Hz, then the PGA for each simulation is calculated as the geometric mean of the horizontal components.

Results
The results show that taking the radiation pattern correction globally improves the results (Fig. 6). However, there are some limitations in the focal areas of the unique Green's function (azimuths -45°and 45°). At an azimuth of 45°, beyond 20 km the correction is no longer taken into account. This is due to the threshold defined (see Sect. 2.5). Indeed, as this is in a nodal area, the unique Green's function coefficients are weak. For an azimuth of -45°, the value of the radiation pattern is even weaker, due to the east dipping orientation of the fault. The correction is then not applied at any distance (Fig. 6).

Empirical Green's functions
Despite the interest generated by the MDF due to its seismogenic potential (see Fig. 1), the seismicity generated by this fault system is low and hence, there are very few events that can be used as EGFs, and they are often limited by their small magnitudes. However, three events were selected here (Fig. 7). The first event was located approximately in the middle of the fault, almost under the city of Manosque. This event occurred on July 8, 2010, at 20:20, with a magnitude estimated to M L = 2.9. This earthquake was the mainshock of a small seismic sequence that lasted for several weeks. The second event was located on the northern part of the fault, near to the town of Villeneuve, and was recorded on September 19, 2012, at 18:56, with estimated M L = 3.5. This earthquake was also followed by a small seismic sequence that lasted for several weeks. This event was, however, too far from the area of interest, and will therefore not be used. The last event occurred on May 9, 2018, at 06:00, with a magnitude estimated as M L = 2.5. This third event was located on the southern part of the fault, near to the CEA Cadarache site. These empirical Green's functions are of relatively low magnitude, consequently the amplitude of low frequency ground motion is relatively small. Therefore, the computed ground motion was high-pass filtered at 1.1 Hz, which corresponds to the minimum limit below which the signal-to-noise ratio is \ 3 (Fig. 8). Note that the EGFs time series have also been visually inspected to ensure the quality of them, including the absence of unexpected pulse waveform dominating the signal (Fig. S1).
Despite the low-frequency limitation, it was possible to estimate their focal mechanisms ( Table 1). The focal mechanism 2010 event was estimated using the FOCmec method (Snoke 2003). The focal mechanisms of the 2012 and 2018 events were estimated using the FMnear method (Delouis 2014). However, the radiation pattern correction , which is only applied up to 3 Hz, will have little impact on these simulations, due to the low frequency limit of the EGFs. Two stations on the CEA Cadarache site are common to the three EGFs, thus limiting the simulations at these two sites: the rock site CA02 with V S30 = 1800 m/s, and the sediment site CA04 with V S30 = 440 m/s ). These V S30 were estimated from in situ measurements involving different approaches: crosshole measurements, surface-wave-based methods, etc. (Garofalo et al. 2016;Ameri et al. 2017a). These sites are equipped with accelerometers: GeosigAC23 for the 2010 event, and GeosigAC73 for the 2018 event. To deconvolute the three EGFs from the influence of their respective sources (Eq. 14), their corner frequencies are also estimated through the Brune (1970) relationship (Eq. 1), using V S = 3500 m/s and Dr = 0.9 MPa, according to the Drouet et al. (2010) inversion of the French Alps data; see the comparison of a theoretical x -2 source spectrum (Brune 1970) with the data spectra (Fig. 8) on the CA02 rock site station and the CA04 sediment site station , sites P2 and P3, respectively). All of the EGF metadata used are summarized in Table 1. The global inversion of the French Alps data realized by Drouet et al. (2010) is intented to determine source terms for each recorded event used, site terms of each station, but also propagation terms such as those describing the geometrical spreading and the anelastic attenuation. Indeed, according to Drouet et al. (2010), V S = 3500 m/s is used for the time correction (Eqs. 15, 16), c = 1.06 for the geometrical spreading correction (Eq. 17), and Q 0 = 336 and a = 0.32 for the anelastic attenuation correction (Eq. 18).

Figure 8
Comparisons of the theoretical x -2 source spectra with the data spectra of the CA02 rock site (a, c) and the CA04 sediment site (b, d). Theoretical x -2 source spectra are computed for Mo and Fc listed in Table 1. The data spectra are also compared with the noise levels (black curves). The low frequency limit, defined by signal-to-noise ratio \ 3, is shown in green. This low frequency is: 1.0 Hz (a), 1.

Geometry of the Fault
The geometry of the fault is built according to Guyonnet-Benaize et al. (2015). The fault is represented as a rectangular area 28 km long by 15 km wide, which is deformed according to the strike and dip variables. The strike progressively varies from 34°to the south, to 13°to the North, and it does not vary with depth ( Fig. 9, left). The fault has a listric character toward the West, with strong dip to the East that decreases to the West. However, this is highly variable in the sense that the strike is not constant, and the dip does not vary progressively according to the depth (Fig. 9, right). This fault is then divided into subfaults the dimensions of which are fixed so that the simulations are valid up to a maximum frequency f kmax = 35 Hz. Thus, according to Eq. (6), the subfault dimensions are fixed to SF dim = 35 m.

Ground Motion Variability
A M 6.0 magnitude event is generated on the southern part of the fault (i.e., close to the CEA Cadarache center). The variability of the ground motion is assessed through generation of a set of 100 simulations where the source and rupture parameters vary for each of them ( Table 2). The corner wavenumber k C in Eq. (3) is fixed according to the theoretical corner frequency f C (Eq. 4), which is therefore directly dependent on the stress drop (Eq. 1). Then, the global ground-motion variability of the simulations is constrained by the stress drop variability. This is constrained according to a lognormal probability distribution (Fig. 10), with a mean of Dr = 0.9 MPa (Drouet et al. 2010), and a standard deviation of s = 0.2923, according to Berge-Thierry et al. (2003) that remains comparable to those of more recent GMPEs. The GMPE sigma at 34 Hz (the frequency associated to the PGA) represents the interevent variability, as it is similar to the intrinsic variability of the data used for regression. The Berge-Thierry et al. (2003) (BT2003) GMPE is considered here, as it is the one that is recommended in the French reference document to assess seismic hazard in the nuclear safety context (Fundamental Safety Rule RFS-2001RFS- -01 2001. The results of this study will be compared with this BT2003 GMPE established in 2003 and with more recent ones where the site classifications are more consistent with the shearwave velocity characteristics of the Cadarache stations studied (see Sect. 4.4).     The stress drop varies between 0.41 MPa and 2.0 MPa (Fig. 10). The rupture velocity V R commonly varies between 0.7*Vs and 0.85*Vs (Heaton 1990). Then, V R is defined according to a uniform distribution between these values. For the 100 cases here, a minimum value of V R of 2450 m/s is obtained, with a maximum value of 2975 m/s. Following Hisada (2001), rupture-time perturbations are also introduced, defined according to a k -2 model (Eq. 8). The characteristic dimensions of the rupture-time variation are randomly fixed between 30% and 70% of L and W (Eq. 9), and the rupturetime variation values are normalized so as not to exceed 10% of the theoretical rupture time (Fig. 11). This value of 10% chosen for the normalization is relatively weak, to avoid the generation of high frequencies in addition to those provided by the k À2 model. The nucleation point is randomly defined for each simulation in the deepest half of the rupture area (Figs. 11, 12), so as to be consistent with Mai et al. (2005). The variation of the static slip, the nucleation position and the perturbation of the rupture velocity, make it possible to generate a set of absolute source functions that are highly variable but that nevertheless remain close to the x À2 model of Brune (1970), as intended by the method (Fig. 11). Note however that the x À2 decay is expected only between f 1 ¼ 1:10 Hz and f max = 20.5 Hz and is supposed to be as x À1 before f 1 .
The rupture dimensions are defined according to the stress drop and the rupture velocity, and in the 100 simulation cases this varies between L Â W ¼ 10360 Â 5600 m and L Â W ¼ 17710 Â 9555 m (see different examples in Figs. 11 and 12). As the dimensions are derived according to V R , its variability has no impact on the corner frequency of the simulated event. However, the variation of the dimensions remains a source of variability, as the ground motion amplitude saturation phenomenon depends on it (Yenier and Atkinson 2014).
Finally, the last source of variability comes from the position of the rupture area on the MDF. The rupture area position is defined so as to be on the southern part of the MDF. First, a variation of the rupture area position varies the minimum distance between the rupture area and the station (from 6.3 to 13.5 km for CA02 station, and 5.5-12.7 km for CA04; see Fig. 12). Moreover, the modification of the rupture area position also varies the focal mechanism (i.e., strike and dip) of the simulation event, as the strike and dip vary with the fault geometry (Fig. 12).

Results
The simulations performed are bandpass filtered between 1.1 Hz and 34 Hz (see synthetic for simulations 8, 12, 36 and 52 in Figure S2). The lowfrequency limit is set according to the signal-to-noise ratio of the EGFs (Fig. 8). As the low frequency part is missing in the simulation, i.e., between the corner frequency (* 0.16 Hz for M 6.0 with Dr = 0.9 MPa and V S = 3500 m/s according to Eq. 1) and the lowfrequency filter limit, we compare the simulated results in terms of the PGA only, with various GMPEs for the horizontal motion (Fig. 13, 14). The peak ground velocity (PGV) is indeed more sensitive to the frequencies around the corner frequency (see Dujardin et al. 2015), and our simulations can therefore not be considered as reliable in terms of the PGV. The PGA of each simulation is calculated as the geometric mean of the horizontal component and represented as a function of the hypocentral, Joyner-Boore and rupture distance (Fig. 13). For comparison, we used several GMPEs. First, we used the Berge-Thierry et al. 2003, (BT2003) since this  GMPE is used within the framework of the French regulation for nuclear facilities. This GMPE has to be use with surface-wave magnitude (Ms). Within the context of Provence seismicity, a specific analysis showed that the Utsu (2002) relationship between Ms and Mw is applicable. In this case, the approximation Ms * M is pertinent around magnitude 6. The BT2003 GMPE uses the hypocentral distance (R HYP ) metric and can account for two site conditions ('rock' for V S30 [ 800 m/s and 'soil' for 300 m/s \ Vs \ 800 m/s. In addition to BT2003, we also compare simulation results with the following more recent GMPEs: Ameri et al. (2017b), Akkar et al. (2014), Bindi et al. (2014), Boore et al. (2014), Cauzzi et al. (2015), Campbell and Bozornia (2014), Abrahamson et al. (2014), Chiou and Youngs (2014) and Laurendeau et al. (2017). These GMPEs will be referred respectively as AM2017, AK2014, BI2014, BO2014, CA2015, CB2014, AB2014, CY2014 and LA2017. All these GMPEs use moment magnitude Mw and will be used here for strike-slip source mechanism. The AM2017, AK2014, BI2014 were derived using a Euro-Mediterranean database. BO2014, CA2015, CB2014, AB2014 and CY2014 were derived using a global database. LA2017 was derived using a subset of the Kiknet database, with a very different approach that will be discussed latter. AK2014, AM2017, BI2014 and BO2014 are computed using the Joyner-Boore distance (R JB ). AB2014, CA2015, CB2014, CY2014 and LA2017 are computed using the distance to rupture plan (R RUP ). The AM2017 is implemented for 4 EC-8 soil classes, other use the real value of V S30 parameter within a given velocity range.  Figure 13 shows the results of simulations for stations CA02 (rock site) and CA04 (soil site) for three GMPEs expressed in 3 different metrics: one in R HYP (BT2003, Fig. 13a), one in R JB (BO2014, Fig. 13b) and one in R RUP (CB2014, Fig. 13c).
First, we note strong variability in terms of distance for a given station (e.g., from R HYP = 9.8 km to R HYP = 21.7 km for the CA02 site, and R HYP = 9.0 km to R HYP = 20.6 km for the CA04 site; Fig. 13a). Such variability is only possible due to the extended simulation code, which allows variations of the size, the position of the rupture area, and the position of the nucleation. However, despite numerous sources of variability, the ground motion variability is in agreement with that of the GMPEs. This is expected since the simulation variability is constrained according to BT2003 variability. We obtain a mean PGA of 0.4820 m/s 2 for CA02, and 3.1465 m/s 2 for CA04, which defines an amplification factor of 6.5 between CA02 and CA04. This amplification is attributed to the strong high frequency site effects on the CA04 station, with amplifications of up to a factor of 10 between CA02 and CA04 for frequencies above 2 Hz . This difference is perhaps greater than what would be observed for real data, as these simulations are only based on high frequencies ([ 1.1 Hz), which correspond to the frequency range affected by the site effects. One can note that the nonlinear response expected for a large event is usually underestimated when simulations are performed with EGFs, even at rock sites.
We note that the BO2014 is a rather good approximation for the CA04 sediment site (Fig. 13b) with a median very close from the mean of simulated PGA values. For BT2003 and CB2014, EFG simulations shows a mean acceleration few tens of percent higher that the GMPE median (Fig. 13a, c), and the GMPEs largely overestimates the ground motion on the CA02 rock site (Fig. 13). Note that this underestimation may also be due to the use of a sparse set of EGFs which is a limitation of the method in this particular MDF application case.
We test more GMPEs on Fig. 14. Figure 14a (resp. 14b) shows results expressed in R JB (resp. R RUP ) for the soil site (CA04). We use the real V S30 CA04 values (440 m/s) for all GMPEs to draw corresponding curves, except for the AM2017 for which we use the 'B' EC8-soil class. Figure 14c (resp. 14d) shows same kind of results for the rock site (CA02). Here, we use the upper V S30 validity limit of GMPEs to draw curves (1200 m/s for AK2014, BI2014, CA2015; 1500 m/s for BO2014, CB2014, AB2014, CY2014), the 'A' EC8-soil class being used for AM2017. On Fig. 14d, we also use the LA2017 as a correction factor applied on another GMPE (here AB2017). Indeed, the LA2017 is derived from a very specific database (Kiknet), it is implemented with quite simple functional form and is built using a quite reduced amount of data. It did not intend to be used 'as it', in particular in an European context. However, it integrates a specific site amplification correction applied station by station before its derivation in order to produce a better estimation for ground motion for hard-rock sites. Here we computed the ratio between the prediction for two site conditions: 440 and 1800 m/s and we apply this correction factor to the AB2014 computed for V S30 = 440 m/s. For the soil site (CA04) we see that most of the GMPEs moderately underestimate the simulation of a few tens of percent. This observation remains consistent with the fact that CA04 is affected by high site Figure 14 Representation of the PGA for the 100 simulations for the sites CA04 (a, b) and CA02 (c, d). Simulations are compared versus the Joyner-Boore distance with AK2014, AM2017, BI2014 and B02014 GMPEs (a, c), and compared versus the rupture distance with AB2014, CA2015, CB2014, CY2015 and AB2014 with LA2017 correction GMPEs (b, d)

2274
A. Dujardin et al. Pure Appl. Geophys. amplification . For the rock site (CA02), on the contrary, all GMPEs strongly overestimate the simulations. The best results (although still significantly overestimating the simulations) are obtained using the LA2017 correction ratio (here applied on AB2014 GMPE). This general ascertainment confirms the need to continue efforts towards the establishment of consistent GMPEs applicable to hard-rock conditions.

Discussion and Conclusions
The EGF simulation code presented here was developed for the specific case of low-to-moderate seismicity areas such as metropolitan France, where the main limit is the number of seismic events that can be used as EGFs. To overcome the lack of recorded data for small earthquake, best use is made of the few available EGFs by shifting them from their original positions to sample the whole fault plane, considering some adjustments. In particular, we implemented a radiation pattern correction, the influence of which is preponderant particularly at the shortest distances. This correction was initially tested and validated in a purely numerical case. The results appear promising, with noticeable improvements at the lowest distances, even if we also note some limitations when approaching the nodal planes of the considered EGFs. However, this test was carried out only for a homogeneous medium and deserves to be tested in more realistic media.
We then use this code to simulation of ground motion from M 6.0 scenario earthquakes on the MDF, where the interest in terms of seismic hazard assessment appears obvious. Indeed the MDF is an active fault with a complex geometry. Its geometry has been extensively studied, which has resulted in a well-constrained three-dimensional geological-geophysical model (Guyonnet-Benaize et al. 2015) that can be used to perform full seismic wave propagation simulations. These latter are of huge interest to account for complex phenomena such as nonlinear interactions that can occur between the geological medium and the buildings under strong seismic motions or to perform full physics-based strong-motion predictions for sites or facilities located very close to the MDF system. Such simulations are however limited at high frequency (above * 1 Hz).
In the present study, instead of using the threedimensional velocity model, we took advantage of the available EGFs, which naturally contain information about propagation and site effects over a broad frequency range. During this application, we encountered the classical limitation of the EGF technique at low frequency: due to the small magnitude of the EGF, the signal-to-noise ratio is too small below 1.1 Hz. This limitation had two consequences. First, we did not compute the PGV values, which are Vol. 177, (2020) Optimization of a Simulation Code Coupling Extended Source (k -2 ) mainly sensitive to frequencies around the corner frequency of the simulated event (fc = 0.16 Hz on average) while our low frequency limit is fmin = 1.1 Hz. However, this limitation has a negligible impact on the PGA (e.g., Dujardin et al. 2015). Second, the proposed correction of the radiation pattern, which is effective up to 3 Hz only, has a very weak impact on the simulated acceleration. This raises the question of how to properly model the low frequency ground motion in the framework of EGF simulations. A first solution is the use of EGF recorded by velocimeters rather than accelerometers. We made a comparison for the 2018 event for a site in the CEA Cadarache (CA01 site; which was not used in the present simulations), which is equipped with both an accelerometer (GeosigAC73) and a velocimeter (Guralp CMG3T) (Fig. 15). The velocimeter has a good signal-to-noise ratio from 0.3 Hz (Fig. 15b), while the accelerometer performs well only above 1.3 Hz (Fig. 15a). Another solution is to complement the EGF simulations with low-frequency three-dimensional numerical simulations (in case there is enough information about the propagation medium) or Green's functions using noise correlation (e.g., Denolle et al. 2013Denolle et al. , 2014. Furthermore, the proposed EGF simulation code accounts for complex seismic rupture processes. First, it is possible to incorporate variable strike and dip angles to model complex fault systems as the MDF. Second, the code models complex rupture propagations and generates populations of rupture scenarios with variable rupture dimensions, positions of the rupture area on the fault, slip and rupture time distributions, and nucleation positions. Note that the variability of the simulations can be mastered by constraining the stress drop distribution, which is an input of the code. This makes it possible to generate a ground motion variability in agreement with that of the accelerometric databases, which is an advantage for seismic hazard assessment studies (e.g., for comparisons between site-specific simulations and GMPEs).
Finally, our results are compared with various GMPEs. GMPEs moderately underestimate simulations at the CA04 sediment site (V S30 = 440 m/s) and strongly overestimate the CA02 hard rock site simulations (V S30 = 1800 m/s). As mentioned before, the recording stations of the Cadarache site benefit from direct in situ measurements that enable them to be assigned very well constrained V S30 . This is a key information since most of the strong-motion records available in databases worldwide (on which the GMPEs are based and obtained) are not yet associated to a well characterized station, with respect to its specific geological and geophysical properties. These latter properties (e.g., geological stratigraphy, density, P wave and S wave-velocities, depth of main impedance contrasts) are directly responsible for the local seismic site responses.
It is also necessary to point out some important aspects about the use of GMPEs (as detailed in Berge-Thierry et al. 2017a, b, concerning the BT2003 with respect to its use within the framework of the French nuclear regulation). Using generic GMPEs to predict ground motion for a specific site where V S30 is far from the class interval boundaries necessarily provides inappropriate results. From this point of view, the results can be anticipated for comparisons between the physics-based ruptures and site-specific predictions of the present study and the GMPEs study. Indeed, with V S30 of 440 m/s, CA04 is clearly in the subset of records that are well represented in the motion databases used to derive GMPEs and their associated soil-site coefficients. On the contrary, with its particularly high shear-wave velocities, CA02 is out of the range of the data that are considered to be representative of ''hard-rock'' responses. An accurate strong-motion prediction for such high V S30 is not possible using generic GMPEs (especially if their functional forms do not directly integrate V S30 as a proxy). This highlight the need of continuing effort of development of alternative approaches as the one proposed by Laurendeau et al. (2017).
These observations raise several points. First, it raises the question of the ground motion on the hardrock site, as there is no clear nor consensual definition of the rock, which is particularly problematic as it is needed to get a reference at least to assess the additional specific site effects coming from: • weathering factors, even on originally healthy rock, that induce lower velocity values within the shallowest layers (first few meters) beneath the surface (see e.g., Hollender et al. 2017). These effects can induce large amplifications of ground motion at high frequency (e.g. [ 10 Hz), even on rock and hard-rock sites; • impedance contrast or soil properties; • finally, geometric configurations and fillings for specific substratum.
The GMPEs are indeed based on the use of real databases where characterization of the so-called rock sites is still poorly, if at all, constrained with very few measurements of the V S30 proxy and/or the frequency resonance of the site effects. The average prediction is therefore naturally higher than what is expected for a true hard-rock site. Secondly, this result raises the question of taking into account the site effects. Although GMPEs are not suitable for site-specific predictions, the improvement of the functional forms, however, allows at least the onedimensional site effects to be taken into account. Indeed, some GMPEs consider V S30 and/or a site frequency of resonance as a proxy of the functional form. Despite this, more complex effects, such as two-dimensional or three-dimensional resonance basin effects (Causse et al. 2009;Dujardin et al. 2016), or local topographic effects due to short wavelength focusing are still not possible. Consideration of such complex effects must therefore be carried out differently, based, for example, on sitetransfer responses.
Finally, the results raise the question of the distance definition used. It appears intuitive to think that R RUP is better suited for the extended fault. This is especially true when distances are small. The results indeed appear to be more consistent when these are represented according to R RUP or alternatively to R JB .