Stochastic Generator of Earthquakes for Mainland France

: Mainland France is characterized by low-to-moderate seismic activity, yet it is known that major earthquakes could strike this territory (e.g., Liguria in 1887 or Basel in 1356). Assessing this French seismic hazard is thus necessary in order to support building codes and to lead prevention actions towards the population. The Probabilistic Seismic Hazard Assessment (PSHA) is the classical approach used to estimate the seismic hazard. One way to apply PSHA is to generate synthetic earthquakes by propagating information from past seismicity and building various seismic scenarios. In this paper, we present an implementation of a stochastic generator of earthquakes and discuss its relevance to mimic the seismicity of low-to-moderate seismic areas. The proposed stochastic generator produces independent events (main shocks) and their correlated seismicity (only aftershocks). Main shocks are simulated ﬁrst in time and magnitude considering all available data in the area, and then localized in space with the use of a probability map and regionalization. Aftershocks are simulated around main shocks by considering both the seismic moment ratio and distribution of the aftershock’s proportion. The generator is tested with mainland France data.


Introduction
Mainland France seismicity is considered as low to moderate due to its high return periods and weak maximal magnitudes. Only a few earthquakes have produced damages in this territory since the introduction of the French insurance compensation system in 1982. The six major earthquakes occurring after 1982 resulted in less than EUR 600 M of insured losses, representing approximately 1.5% of the total amount compensated with this compensation scheme, https://catastrophes-naturelles.ccr.fr/ (accessed on 20 December 2021). Nevertheless, major earthquakes could strike this territory (e.g., Liguria in 1887, Basel in 1356) and create financial losses and casualties. For example, CCR (Caisse Centrale de Réassurance, French insurance company) and BRGM (Bureau de Recherches Géologiques et Minières, French geological survey) have quantified the probable insurance losses associated with historical earthquakes if they occurred in the present day: It shows that estimating the seismic hazard and risk is necessary even in low-tomoderate seismic areas, since losses can be significant.
The Probabilistic Seismic Hazard Assessment (PSHA) [5] is a group of methodologies that estimates seismic hazard in a probabilistic way using seismic, tectonic and geologic data, as well as their uncertainties. It leads to a measure of probabilities of exceedance for different hazard levels (e.g., Peak Ground Acceleration) in a given site. The spatiotemporal analysis of past seismicity recorded in an area leads to distributions of the earthquake's occurrence and location that are the cornerstone of the PSHA. Traditionally, these distributions can be used in an integrative way [5] or as an input to constrain a stochastic generator of earthquakes [6]. The main advantage of the latter is to explicitly consider some synthetic seismic events that allow for the computation of the contributions of these earthquakes to the hazard evaluated at a given site. This aspect is useful in an insurance context and also makes it possible to link earthquakes to others hazards (e.g., [7] for tsunamis, [8] for liquefaction, [9] for landslides and so on). Regardless, this study only focuses on earthquakes.
The PSHA process takes into account the spatio-temporal behavior of seismicity. To carry this out, three principal methods can be involved: • Homogeneous seismotectonic zones, hereafter called the zoning method [5]. Seismicity and tectonics are considered homogeneous within each zone. Seismic temporality is defined by computing a specific Frequency-Magnitude Distribution (FMD) in each zone. These zones are produced from various and complex data, such as seismic, geodesic and geologic ones, e.g., [10] (for France application); • The smoothing method [11,12]. FMD is computed smoothly in space based on the past seismicity; • The active faults method, e.g., [13] (for recent California application). Each known active fault is analyzed in order to associate it with its own FMD.
Each of these methods have drawbacks and advantages. These methods are classically used in parallel in order to take into account the advantages of each. Their application depends on the user's wishes and mainly on the studied region.
Mainland France is a low-to-moderate seismic area since it is localized on an intraplate domain, which is also called a stable region, i.e., far from tectonic plate boundaries.
The spatial segregation of the zoning method induces a drastic reduction in data. In fact, available data in each region are far smaller than those observed in the whole studied area. This finding is even more critical in low-to-moderate seismic areas, where seismic data are scarce. Using this method calls for specific care regarding data completeness and representativeness.
The smoothing method provides a good description of the seismicity but is not complete. In intraplate regions, where seismic frequencies are low, longer observation periods are needed. This is particularly true for extreme events (highest magnitudes) that may not have been seen during the period of observations. In fact, according to some authors, e.g., [14,15], stable regions present a seismicity that is more diffuse and homogeneous than past observations. Finally, in intraplate domains, only few or even zero active faults are known. Moreover, a seismic cycle assumption could be questionable in stable regions [15]. According to these authors, intraplate seismicity is not driven by repeating tectonic loading and release energy but by transient crust deformations [15] caused by other origins than tectonics; see [16] (for a review in a French context). For this reason, the application of the active fault method is disputed. For example, for an application in the southeastern part of France [17], Martin et al. used this method as a complement to others, but only for the Provence region, where active faults are better known, e.g., [18,19].
In this paper, we propose a new stochastic generator to simulate the seismic catalogues of earthquakes that considers the specificities of the intraplate region (e.g., little data, various and complex origins of seismicity). We start by presenting the stochastic generator that simulates both main shocks (main and independent events) and aftershocks (correlated events). Then, we apply this methodology to mainland France.

Methodology
The proposed generator first simulates main shocks and then their associated aftershocks. A global scheme of the method is illustrated in Figure 1. In this scheme, the method is represented by the blue blocks, whereas inputs are marked as green blocks. Data used to produce these inputs are referred to in orange. Since this paper only focuses on methodology, the user is free to use their own data.

Main Shock Generation
Only few data are available in low-to-moderate seismic regions. Knowing this, instead of computing a Frequency-Magnitude Distribution (FMD) per zone or a smoothing FMD, we established the FMD of main shocks recorded on the whole territory. This has the advantage of maximizing the number of data and evaluating the probability of occurrence of an earthquake in France without any prior assumption. This follows recent studies, e.g., [14,15], that stipulate that stable region's seismicity is more diffuse and homogeneous than past observations suggest. Main shocks are thus generated through time (year) and magnitude for the whole territory. Then, the position of the main shock is moved according to its magnitude and a fault probability map.

Main Shock Generation in Time and Magnitude
Synthetic main shocks are first generated through time and magnitude. The time step chosen is the year, while the magnitude step is set to 0.1. We first set the number of years N y through which the synthetic catalogue will span. Then, the generator calls for a stochastic Frequency-Magnitude Distribution (FMD) of main shocks, which corresponds to a probabilistic distribution of annual number of main shocks for each magnitude step. The annual density λ y (M) of main shocks is defined year after year y for each magnitude step M according to this stochastic FMD of main shocks.
Main shocks are simulated yearly following a homogeneous Poisson point process of annual density λ y (M): with n y (M) being the number of main shocks of magnitude M ∈ [M min , M max ] generated at the year y ∈ [1, N y ]. The annual number of main shocks of magnitude equal to M at the year y (λ y (M) Equation (1)) is drawn in the stochastic FMD of main shocks. However, these FMDs are defined for main shocks with magnitude greater than or equal to M (λ y (≥M)). The annual density obtained from this kind of FMD is not λ y (M) but λ y (≥M). Thus, λ y (M) is defined as: In some cases, λ y (≥M + dM) can be higher than λ y (≥M), and so λ y (M) is set to 0. This step is referred as 1A in Figure 1.

Main Shock Generation in Space
Every main shock simulated in terms of year and magnitude is then localized somewhere in the whole of mainland France.
A recent study [16] stipulates that the structural inheritance-that is, crust or mantle weakening from past tectonic-can play an important role in deformation localization in intraplate context. According to previous works and new observation, these authors estimate that 55-95% of the stable region's seismicity occurs near such weaknesses. Based on this, we suppose here that such geological objects, that give testimony to hundreds of millions years of seismicity, give more exhaustive information than past seismicity in predicting possible location of future earthquakes. In this paper, for a first approximation, we propose guiding the location of generated main shocks by using a fault probability map.
However, even in active regions around the world, fault activity could be difficult to address, e.g., [20] (for the Canterbury earthquake in New Zealand). Recent studies in China [21] or earthquakes in France [22] and Australia [23] serve as reminders for how difficult it is to define fault activity in stable regions. Since, in these regions, faults can be inactive over long periods of time [14,21], we decided to consider all faults, regardless of their presumed importance and activity.
The map of faults is converted into a density map by computing the average length of fault lines per unit area on a Cartesian grid. After that, the density values are normalized in order to obtain a sum equal to 1, which leads to a probability map. In that sense, each cell of the grid is characterized by a probability of hosting faults. Grid resolution is set by the user.
Looking at past seismicity, one would notice that earthquakes with the highest magnitudes are concentrated in the most active regions, even in intraplate domains. For example, these regions in France are the Alps, the Pyrenees and the Rhine basin. We decide to apply a magnitude threshold to drive highest magnitudes' earthquakes in areas where they have most chance to happen. To achieve this, we segregate the fault probability map by regionalization. Each region is characterized by a maximal magnitude M max region . This magnitude corresponds to that of the largest earthquakes that this region can host.
Coupling regionalization with probability map leads to a set of regions defined by their own probability map and maximal magnitude M max region . A main shock of magnitude M generated at year y can only be localized in eligible regions, i.e., regions where M ≤ M max region . Once the set of eligible regions is defined, the location of this main shock is drawn through the cumulative distribution function of the corresponding probability map ( Figure 2). This procedure is a heterogeneous Poisson process whose density depends on the magnitude of the synthetic main shock to be located.
The Cartesian grid on which fault probability map is computed is used to select a cell, whose resolution depends on the user, where the main shocks' coordinates (the black point in Figure 2) are uniformly drawn. Figure 2 summarizes the localization process. Owing to this approach, spatial distribution of generated main shocks is limited by magnitude, with respect to the fault probability map that characterizes seismicity.
This step is referred as 2A in Figure 1.

Rupture Plane's Parameters
Seismic hazard produced by an earthquake at a given site is computed by groundmotion prediction equations. These equations accept the definition of earthquakes as points. However, the geometry of rupture plane associated to an earthquake could play a significant role in seismic hazard computation. Knowing this, a rupture plane needs to be described for each generated earthquake.
In this study, a rupture plane is considered as a plane at a particular depth and is parametrized by a length (L) and three orientation angles: (i) azimuth (clockwise angle from north between 0 • and 360 • ), (ii) dip (angle between 0 • and 90 • from the horizontal) and (iii) rake (angle between −180 • and 180 • ) that define the components of fault movement (normal, reverse, strike-slip).
L is evaluated from moment magnitude M w of the earthquake using the following relation, e.g., [24,25]: with L being the rupture plane's length in km and l 1 and l 2 being two constants. For the other parameters, they are defined per region, where intervals of values are set according to data. Classically, these kinds of data are used: • Instrumental seismic catalogue for depth; • Focal mechanism catalogues for azimuth, dip and movement; • Fault map for azimuth.
See Appendix A for an example of definition of these intervals for an application in mainland France.
Values for each parameter are randomly drawn in a uniform distribution when a main shock in generated in a given region.
This step is referred as 3A in Figure 1.

Aftershock Generation
Aftershocks could play a significant role in the seismic risk, e.g., [26], since their occurrence can destroy a building damaged by the main shock. Generating aftershocks is thus necessary in the context of seismic risk. Once again, the low amount of data makes seismic sequences difficult to observe in intraplate domains. Thus, applying the famous Epidemic Type Aftershocks-Sequence (ETAS [27]) model, which is a marked point process fond of parameters, is difficult.
Computing Frequency-Magnitude Distribution (FMD) of main shocks in order to generate main shocks through time and magnitude calls for differentiation between main shocks and aftershocks in the data analyzed. This leads to a Proportion-Magnitude Distribution (PMD) of main shocks that depicts the proportion of main shock in the studied catalogue in function of magnitude. Here, we propose generating aftershocks complemen-tarily to main shocks by following the same PMD used to produce the main shocks. Once every main shock is generated, aftershocks are produced and related to their own main shock. This aftershock's production is carried out through three steps described below and represented in Figure 1 (1-3B).

Number of Aftershocks to be Produced
The first step consists of defining the number NbAs(≥M) of aftershocks with magnitude greater than or equal to M to be produced. This number is obtained for each magnitude step according to the number of main shocks NbMs(≥M) already generated and the proportion of main shocks PropMs as follows: An example of PMD of main shocks, as well as the number of main shocks with magnitudes greater than of equal to 4 generated over 100,000 years, are shown in Figure 3. The number of aftershocks M ≥ 4 to be produced, computed by Equation (4), is also visible in Figure 3. Thus, aftershocks are firstly associated with their magnitude.   (4)).

Aftershock-Main-Shock Relation
According to Båth law [28,29] (Equation (5)), main shock is the event of the sequence with the highest magnitude. Moreover, its magnitude M ms is limited from below according to the magnitudes M as of its aftershocks as follows: where ∆M is a constant equal to 1.2 for crustal seismicity (<50 km depth, in [28]).
In this paper, we use this convenient law to find the eligible main shocks for each aftershock. In fact, owing to the first step, the number of aftershocks to be produced is known for each magnitude step ( Figure 3). Applying Equation (5) for each aftershock allows us to define the minimal magnitude M ms min of its related main shock. Its linked main shock is then randomly selected among all of the main shocks with magnitude greater than or equal to M ms min . Aftershock's year of occurrence is set to the same as its main shock. However, one needs to keep in mind that Equation (5) is an empirical model and thus not a general truth. For example, the M6.5 Amatrice (Italy) earthquake provoked an aftershock of magnitude 6.1 three months later. In order to take into account ∆M variability in the aftershock-main-shock association, the method of seismic moment ratio [30] is used in this paper. It computes this term: according to the R ratio equal to: where ∑ M 0 as and ∑ M 0 f s are the seismic moment, equivalent to energy, released by all of the aftershocks and foreshocks, respectively, and M 0 ms is the seismic moment released by their main shock. According to [30], R is equal to 0.05. A ∆M distribution is obtained by sampling R in a Gaussian law with mean equal to 0.05 and a standard deviation of 0.0125 [30]. Such a distribution gives an average ∆M equal to 0.87 while 90% of ∆M values between 0.77 (Q5) and 1.02 (Q95). Drawing a ∆M in this distribution makes it possible to estimate M ms min , the minimal magnitude of its main shock, from Equation (5).

Aftershock's Location and Rupture Plane's Parameters Definition
We suppose that aftershocks are localized on the same fault as the main shock. Aftershocks' epicenter (x as , y as ) are thus placed at approximately 0.75 L from main shocks' epicenter (x, y) in the direction of its azimuth (azi) to within 10 • : where L (km) is the rupture plane length associated with the main shock (Equation (3)). Depth z as (km), azimuth azi as ( • ) and dip dip as ( • ) of aftershock's rupture plane are sampled into normal laws. These laws have a mean that is equal to the value of the main shock (z, azi et dip) and a standard deviation depending on the parameter: Movement of aftershock's rupture plane (normal, reverse, strike-slip or unknown) is equal to the main shocks' rupture plane.

Regionalization
Even in intraplate domains, the spatial distribution of seismicity can seem heterogeneous. For example, in France, the Pyrenees and Alps are linked to a collision of European and African plates, which accentuates the tectonic features and seismicity in these areas. The SHARE project [31] has proposed segregating European territories according to tectonic features by differentiating pure Stable Continental Regions (SCR) and Oceanic Crusts (OC) from shallow active regions. This tectonic distinction is visible in Figure 4 ( Figure 2 in [32]) for France. This regionalization is used in this paper in order to limit the spatial distribution of magnitudes (see Section 2.1.2 and Figure 2).

Seismic Catalogues Historical Catalogue
The historical French catalogue FCAT-17 [33] is composed of 4250 earthquakes felt and reported by the French population from 463 to 1964. Despite being associated with large uncertainties, such a long catalogue is useful to define the maximal magnitude of the study (M max ) and per region (M max region ). In this study, we decide to select these maximal magnitudes as the maximum recorded magnitudes in the FCAT-17 catalogue, plus its uncertainty. The strongest earthquake of this catalogue occurred in Liguria in 1887, with a moment magnitude of 6.7 ± 0.6. Thus, the maximal magnitude M max is set to 6.7 + 0.6 = 7.3. Since this earthquake is localized in region 4 (compressional Alps, Figure 4), this region is characterized by a M max region equal to 7.3. Table 1 summarizes the maximal magnitudes observed in each SHARE region. Consistent with the recent study [34], we consider that a seismic event with a magnitude of up to 5.5 can occur everywhere in France. As a result, the M max region is raised to 5.5 if the estimation using the catalogue leads to a smaller value, which is the case for the OC region (n • 2, Table 1). Maximal magnitudes are constant and, thus, their uncertainties are not explored in this paper.

Instrumental Catalogue
This study uses both the instrumental French catalogue SIHex [35] and the RéNaSS catalogues, https://api.franceseisme.fr/fr/search (accessed on 20 December 2021). The former is composed of 37,408 earthquakes recorded between 1965 and 2009, whereas the latter regroups 25,042 seismic events indexed from 2010 to 2020. In total, the merged catalogue brings together 62,450 earthquakes that occurred during the 1965-2020 period.
Such a catalogue isn't complete since the seismometers' number and resolution have improved over time. Also, the catalogue needs to be processed in order to reveal representative and meaningful information. The classical approach to make this catalogue exhaustive is to consider cut-off magnitudes (M c ) and years (Y c ). The catalogue can be considered exhaustive on a particular territory for earthquakes with a magnitude greater than or equal to M c only during the period Y c -2020. For the French territory, we determine M c and Y c thanks to the cumulative visual method [36,37]. The couples of M c and Y c that we found are visible in Table 2.

Faults
For a first approximation, we propose building a spatial probability map ( Figure 2) by using the CHARM database, https://infoterre.brgm.fr/page/telechargement-cartesgeologiques (accessed on 20 December 2021), produced by BRGM (Figure 6a). This database contains fault traces, i.e., intersection lines between faults and the Earth surface mapped on geological maps. The CHARM database is composed of 10,576 terrestrial fault lines covering the whole of mainland France and its vicinity (Figure 6a). As already stated in Section 2.1.2, we decide to consider all faults, regardless of their presumed importance and activity. This map of fault traces is converted into a density map by computing the average length of fault lines per unit area on a Cartesian grid of 5 km. The density map is illustrated in Figure 6b. a b Spatial distributions of past earthquakes ( Figure 5) and fault traces ( Figure 6) show the same trend. However, one can remark some differences between these two maps. For example, some moderate earthquakes have been observed in the north of Bordeaux, whereas no fault trace is visible nearby. These differences are explained by: • Unknown, e.g., [20], and/or new faults; • The fact that the CHARM database only contains terrestrial fault traces, which means that covered, deep and/or bathymetric faults are missing; • Very old faults that have probably been inactive for a long period of time (for example, Armorican Massif).
In order to allow for the possibility of localizing main shocks anywhere [14,15], thus avoiding specific cases, such as the one observed near Bordeaux, regions where no fault lines have been documented are associated with a minimal density value instead of 0. In this study, in the first instance, this minimal value is arbitrarily set to 1% of the maximal density value. In that sense, we consider that the least represented regions have a density that is 100 times lower than in most represented regions. Concretely, according to Figure 6b, every 5 × 5 km cell characterized by a density value lower than 0.34 100 = 0.0034 is associated with this minimal value (0.0034). The cells concerned by this density modification represent around 15% of mainland France. It mainly concerns Parisian and Aquitaine basins (SCR), as well as the Mediterranean Sea (oceanic crust). This map is then transformed to a probability map due to normalization. The final regionalized fault probability map is used as an input of the generator of earthquakes for mainland France ("Regionalized density map of faults depending on magnitude", Figure 1).
Faults (Figure 6a) can also be employed to define the probability density functions (pdf) of azimuth in each region. Please refer to Appendix A for an application to mainland France.

Frequency-Magnitude Distribution
The frequency-magnitude distribution (FMD) describes the temporal occurrence of earthquakes. It consists of computing the annual number of earthquakes as a function of the magnitude. The classical way to compute FMD is to apply the GR law [38] on observed data: where M is the earthquake's magnitude and N(≥M) is the annual number of earthquakes with a magnitude greater than or equal to the M observed in the seismic catalogue. In this study, M stands for the moment magnitude (M w ), and it is limited by a minimal magnitude M min (which, here, is equal to 2). The seismic catalogue is composed of two types of earthquakes: main shocks and their correlated seismicity (fore shocks or aftershocks). Since they are characterized by different spatial and temporal behaviors, seismologists usually analyze them separately. A declustering algorithm is classically used for that purpose, e.g., [39]. It applies spatiotemporal windows around earthquakes in a catalogue in order to find clusters. In each cluster, the earthquake with the maximal magnitude is considered as the main shock, and the others are considered as correlated events. Thus, it is possible to segregate independent events from dependent ones by conserving only main shocks and alone earthquakes (clusters with only one earthquake). Then, the population of main shocks can be studied separately through a Poisson distribution in time, space and magnitude. In this study, the G85 declustering algorithm [40] is chosen. The Proportion-Magnitude Distribution (PMD) of main shocks obtained by applying this algorithm on French instrumental seismicity is represented by blue empty dots Figure 7.  (10)) on data has been applied from M2 to M5 and has been extrapolated with the use of the truncated GR law (Equation (11)). Figure 7 show the annual number of main shocks observed in the exhaustive and declustered catalogue as a function of their magnitude. Applying Equation (10) to these data requires a straight line distribution. Such a distribution is observed between magnitudes 2 and 5 ( Figure 7). Thus, we apply Equation (10) between M2 and M5 and obtain a = 4.41 and b = −1.12.

Orange empty dots in
From M5.1, the annual numbers of main shocks fall and do not respect the GR law (Equation (10)). This discrepancy for large magnitudes is classically attributed to a noncompleteness of data for these magnitudes. Assuming that, an extrapolation of the GR law that is estimated to be between magnitudes 2 and 5 is carried out for M ≥ 5.1.
As the total energy released by earthquakes is finite, some deviation from the GR straight line is required for largest magnitudes in order to avoid infinite distribution. A truncation is generally applied on the GR law [41]: with β = b × log(10). Applying Equation (11) with a and b estimated before (Equation (10)) gives the FMD of main shocks modeled from the exhaustive and declustered catalogue. This FMD is visible in Figure 7 and is characterized by a mean slope (b value) that is equal to 1.12 and an asymptotic fall of frequencies for the maximal magnitude M max 7.3. The mean return periods (inversely equal to annual frequencies) of main shocks are summarized for six magnitude steps in Table 3. These return periods describe the average waiting time between main shocks of magnitudes greater than or equal to M. As expected, the mean return periods increase as the magnitude grows. Return periods 2.4 ± 0.1 d 31.9 ± 0.3 d 1.14 ± 0.05 y 15 ± 1 y 200 ± 24 y 4700 ± 700 y

Consideration of Magnitudes' Uncertainties through a Monte Carlo Scheme
Taking into account magnitudes' uncertainties in the calculation of FMD through a Monte Carlo scheme allows us to obtain a stochastic set of PMD and FMD that corresponds to the global representation of main shocks' proportion and occurrence (Figure 7). This Monte Carlo process consists of producing 1000 initial catalogues of 15,567 earthquakes with M ≥ 2. Only the magnitude differs from one catalogue to another, since the magnitude of each earthquake is drawn in a normal law N(µ, σ), where µ is its magnitude and σ its magnitude's uncertainty.
In Figure 7, each magnitude step is characterized by a set of annual frequencies of main shocks. Thus, for each magnitude step, the annual frequency of main shocks can be described by a probability density function (pdf). These pdfs are visible in Figure 8 for six magnitude steps. For a given magnitude M, the annual frequency of main shocks of magnitudes greater than or equal to M is thus defined as a pdf. The pdf values are normalized by their sum. In this way, these normalized values are dimensionless and their sum converges to 1, which makes them probabilities.
These pdfs of annual frequencies as a function of magnitude are inputs of the generator of earthquakes for mainland France ("Stochastic FMD of main shocks"), whereas the mean PMD of main shocks constitutes another input ("PMD of main shocks", Figure 1).

Results
This section presents 100,000 years of synthetic seismicity generated by using SHARE regionalization (Figure 4). Generated earthquakes have magnitudes greater than or equal to 4.

Main Shock Generation Cumulative Seismic Moment Distribution
Cumulative seismic moments produced by both synthetic and instrumental catalogues are compared Figure 9. The seismic moment M 0 (N.m) released by an earthquake is calculated for each real or synthetic main shock directly from its moment magnitude M w following [42]: According to Table 2, the completeness period associated with instrumental M ≥ 4 earthquakes is 1965-2020, which represents 56 years. Comparisons are thus made over this completeness period. Since 100,000 years have been generated, we divided the synthetic seismicity into a set of 1785 sub-catalogues of 56 years. Each of these sub-catalogues are comparable to the complete instrumental catalogue.
The final seismic moment observed in the instrumental catalogue is smaller than the final mean seismic moment generated. This slight difference is due to the fact that the return periods of M ≥ 4 main shocks given as an input are equal to 1.14 years (Table 3) and 1.24 years, according to the initial catalogue. Thus, the initial catalogue (red line Figure 9) is composed of 45 M ≥ 4 main shocks over 56 years, whereas the sub-catalogues of 56 years are on average composed of 49 M ≥ 4 main shocks (black line Figure 9). Thus, the generator produces a few more M ≥ 4 earthquakes than we observe in the initial catalogue.  Table 2). The synthetic catalogue used is 100,000 years long.
The mean distribution of generated seismic moments is, by definition, smooth, and could differ from the one observed in the instrumental catalogue, but trends are similar. Distributions visible in the synthetic sub-catalogues show variations, as well as the instrumental distribution.

Frequency-Magnitude Distribution
The Frequency-Magnitude Distribution (FMD) of generated main shocks can be observed in the 100,000 years of generated seismicity. This FMD of main shocks is compared Figure 10 with the mean one given as an input to the generator (Figure 7). The generator manages to reproduce this FMD of main shocks, especially for magnitudes smaller than 6.5. For higher magnitudes, slight variations are observed between the generator's FMD and the data's FMD. This part needs to be improved in order to avoid these slight variations.
The FMD of generated main shocks for the six SHARE regions (Figure 4) are also shown in Figure 10. Table 4 presents the b-values computed by regressing the GR law (Equation (10)) on these FMDs of main shocks, from magnitudes 4 to 5 and between M5 and M6.

Rhine Basin Alps
Extension Data Generator France Figure 10. Mean frequency-magnitude distributions of main shocks given as input (black) and observed in a synthetic catalogue of 100,000 years for the whole of mainland France (grey) and for the SHARE regions ( Figure 4).
Aside from the Rhine basin region (n • 5 Table 4), every region shares practically the same slope of the FMD of generated main shocks before magnitude 5 ( Figure 10): from 1.09 to 1.15. These slopes are close to the slope of the FMD of main shocks given as an input: 1.12. This is expected, since we only use this FMD as an input. Homogeneity is thus observed in the results. Table 5 lists b-values calculated on the data in each region. Please refer to Appendix B for more details on how these values have been computed. Regardless of what is obtained from the observed (Table 5) or synthetic (b(M4-5), Table 4) seismicities, the b-values of a given region keep the same order with respect to the b-value used as the input. In fact, regions 1, 5 and 6 are characterized by both b-values being greater than or equal to the reference b-value (1.12). In parallel, regions 2, 3 and 4 are characterized by both b-values being lower than the reference b-value. Thus, variations of synthetic b-values around the initial value are consistent with b-values obtained from data. Table 5. b-values computed by regressing GR law (Equation (10)) on FMD of main shocks observed in instrumental seismicity for each SHARE region. Region numbers refer to Figure 4. See Appendix B for more details. Moreover, for M ≥ 5.1, slopes of FMD diverge from one region to another. Between magnitudes 5 and 6, the FMD of stable oceanic crust and the Pyrenees are, respectively, characterized by the highest and the lowest slopes: 1.98 and 0.98 (Table 4). The higher the slope, the lower the seismic activity. Thus, according to the generator, stable continental and oceanic regions are the lowest seismic areas, and the most active regions (Alps, Pyrenees and Rhine basin) are the highest seismic areas. These results are in line with our knowledge of mainland France seismicity. The use of an input as a fault probability map and only one FMD of main shocks thus leads to consistent results.

Region
From the magnitude 5.5, the limitation of large earthquakes' spatial distribution through the regionalization seems to be efficient. In fact, when the magnitude is too high to appear in a given region, it moves to another one, which concentrates the generated main shocks in regions where high magnitudes are allowed. Thus, the limitation of large earthquakes' spatial distribution has an impact on the change in b-values, which is noted in Table 4. Moreover, a relatively pronounced asymptotic fall of frequency is noticed at the region of maximal magnitude, except for the Rhine basin region (Figure 10). This observation is also in line with our knowledge of general seismicity: without using the truncated GR law (Equation (11)), the generator seems to reproduce this asymptotic fall at the largest magnitudes. Finally, the sum of each region FMD is equal to the FMD of the whole of mainland France. This behavior is due to the use of only one FMD as an input, and brings coherence to the results. Figure 11 illustrates the number of main shocks and aftershocks per unit of magnitude generated over 100,000 years. The higher the magnitude, the higher the proportion of main shocks. This is consistent with the literature and with the Proportion-Magnitude Distribution (PMD) of main shocks given as an input of the generator (Figure 3), which is identical to the one observed in the synthetic seismicity. Thus, the generator of aftershocks is working well, since it manages to reproduce this PDM. Figure 12 shows the number of aftershocks per main shock as a function of the magnitude of the main shock. This number is constant and equal to 1.6 on average for magnitudes of main shocks between 5.2 and 6.9. Number of generated earthquakes 10 4 Main shocks Aftershocks Figure 11. Number of main shocks and aftershocks generated over 100,000 years as function of magnitude.

On the Use of Fault Probability Map
The use of the spatial probability of faults to guide the location of generated earthquakes is a strong choice. As already stated, it comes from both various seismic origins in stable regions and the fact that existing faults are more likely to localize further seismic events [16]. The lack of completeness and clear definition of active faults have led us to consider all faults, regardless of their alleged importance and activity. This choice is questionable but consistent with the fact that, in stable regions, faults can be inactive over long periods of time [14,21]. Results described in this paper are encouraging, since they are consistent with data given as an input and with our knowledge of French seismicity. However, the method presented in this article corresponds to a first step in a new way to generate stable seismicity, and some improvements must be achieved.
First of all, we have seen that the spatial distribution of past seismicity in mainland France is heterogeneous and that the strongest occurring earthquakes are concentrated in specific regions (Alps, Pyrenees and Rhine basin). To consider this statistical observation, we opt for a spatial limitation of magnitudes by applying regionalization. Still, we can imagine using other methodologies to generate seismicity in the most active French regions. For example, a methodology that is closer to what is conventionally carried out with FMD computed by region/zone or smoothly, e.g., [17,34,43] (for recent mainland France applications). Another possible methodology could be to compute a FMD of main shocks in the Alps, the Pyrenees and the Rhine basin and to use a spatial probability map of past seismicity. One could also use a spatial probability of faults weighted by stress rate or others tectonic and geologic information in order to address the seismic potential of faults.
Furthermore, one can imagine weighting the probability map used in this study with our knowledge of intraplate domains. For example, various origins of stable seismicity in the long term could be considered: topography potential energy, erosion and glacial isostatic adjustment since the last glaciation and so on [16]. Transient processes (e.g., fluid pore pressure increase and hydrological or sedimentary load change [15]) leading to stable seismicity could also be taken into account.
Finally, the use of a fault map calls for exhaustiveness of these faults. However, this is not the case, since the map used is only composed of terrestrial fault traces, excluding covered, deep and/or bathymetric faults. Attempting to complete this database should be one of the next steps of our work.

On the Definition of Maximal Magnitudes
As already stated, 1500 years of historic data (FCAT-17 catalogue) cannot be representative for stable seismicity. That is why we choose to use a fault map instead of an earthquake map in order to drive the synthetic earthquake's location. Thus, our definition of maximal magnitudes is paradoxical, since it is based on historical seismicity.
Analyzing analogical regions in Europe or all around the world through a Bayesian approach, e.g., [44,45], can make the definition of maximal magnitudes more robust and not constant. This approach also has the advantage of defining maximal magnitude not as a constant but across a range of values, which is interesting when investigating maximal magnitudes' uncertainties.

On the Aftershock Production
Epidemic Type Aftershocks Sequence (ETAS) models are marked point processes [27] that produce main shocks and associated aftershocks. For that purpose, they use various laws, such as the aftershocks production law ( [46] from [47]): where K is the number of aftershocks produced by a main shock of magnitude M ms , M c is the minimal magnitude of interest and k and α are two constants. Thus, the number of aftershocks produced by a main shock increases with M ms . In our results (Figure 12), the number of aftershocks per main shock seems to be independent to the magnitude of main shocks. This observation is mainly due to the fact that the Proportion-Magnitude Distribution (PMD) of main shocks used as an input in this study ( Figure 3) is derived from the instrumental catalogue. However, as already stated, this catalogue is far from exhaustive for large magnitudes (M > 5). Although it is not a problem for the FMD of main shocks, since it is produced by the extrapolation of the exhaustive part (see Section 3.2), it is a problem for PMD, which is not extrapolated and is thus not exhaustive. Thus, within the whole M ≥ 4 generated earthquakes, aftershocks represent only 6% (Figure 3), which is low, as Figure 11 illustrates.
A PMD of main shocks computed thanks to the G85 declustering algorithm [40] from the exhaustive historical catalogue (FCAT-17) should be more representative of the true distribution. According to this catalogue, aftershocks represent 15% of the whole M ≥ 4 seismicity. Applying the PMD of main shocks estimated from the exhaustive FCAT-17 catalogue gives results that are shown in Figure 13.
These results seem more realistic, since they are based on more exhaustive data and since the number of aftershocks per main shock follows the Equation (13). Moreover, our method to generate aftershocks is non-parametric in contrast to this equation, which needs to set two parameters (k and α). This is an advantage in the context of low-to-moderate seismicity, where objective parametrization from data is limited since data are sparse.
Main shocks are generated according to the proportion of main shocks p observed in the data, whereas the number of aftershocks produced depends on the complementary of p (1 − p). However, results in Figure 13 have been obtained by using instrumental data to generate main shocks and historical data to produce aftershocks. Thus, two PMDs have been used, and so the consistency between p and 1 − p no longer holds. A solution could be to produce both the PMD and FMD of main shocks with historical data. Nevertheless, this catalogue is characterized by large uncertainties (magnitude, time of occurrence and space) and is known to overestimate magnitudes compared to the instrumental catalogue (e.g., Figure 2c in [16]). Its use must be carried out with care.

Conclusions
In this paper, a new generator of earthquakes is proposed and applied to mainland France. Applying a stochastic generator of earthquakes to a French context is not new (see, for example, Ref. [48] for a Pyrenean application). However, the method of generating synthetic seismicity is new. Classically, two methods exist to compute a Frequency-Magnitude Distribution (FMD) of main shocks: (i) smoothly through a kernel approach [11,12] or (ii) discretely by using a zoning approach [5]. This allows us to analyze spatial and temporal behaviors of seismicity at the same time, but calls for data number reduction. In intraplate domains, i.e., far from tectonic plate boundaries, data are sparse. Thus, contrary to these classical approaches, the proposed generator uses a FMD of main shocks at a national scale in order to generate main shocks in time and magnitude in order to maximize the number of data available. These main shocks are then spatially distributed through the use of a probability map and regionalization. The former is used to guide the location of main shocks, whereas the latter allows us to limit the distribution of large earthquakes in space. Intraplate seismicity seems to be more uniformly positioned than in active regions [14,15], structural inheritance "can play a strong role in deformation localization" [16] and stable faults' activity is difficult to define [21][22][23]. For these reasons, faults, regardless of whether they are supposed as active or not, are used to produce the probability map.
Aftershocks are then produced by using the Båth law [28,29], the seismic moment ratio [30] and the Proportion-Magnitude Distribution (PMD) of main shocks. This approach defines aftershocks and main shocks differently from the famous ETAS models. Nevertheless, it remains consistent, since the numbers of main shocks and aftershocks are complementary and depend on data. Moreover, unlike ETAS models, this approach has the advantage of being non-parametric.
Temporal and energetic behaviors of generated main shocks are in line with inputs (FMD of main shocks and regionalization) and our knowledge of mainland France seismicity.
However, some improvements can be achieved, such as completing the fault map; in particular, by bathymetric faults, and better describing maximal magnitudes in each region. Moreover, we have seen that using instrumental seismicity to produce the DPM of main shocks is not representative enough. Using the FCAT-17 historical French catalogue could be a solution, but it needs to be carried out with care due to its large uncertainties and its overestimation of magnitudes. Finally, the method proposed in this paper should be applied to other stable continental regions, such as Northwestern European countries, Australia and so on in order to test its effectiveness.
Author Contributions: Conceptualization, P.T., C.G. and F.B.; methodology, C.G., P.T. and F.B.; data curation, C.G. and P.T.; writing-original draft preparation, C.G.; writing-review and editing, C.G., F.B. and P.T. All authors have read and agreed to the published version of the manuscript.
Funding: This work was performed in the frame of the RING and EXTRA& CO/ICEEL projects at Université de Lorraine. It has been funded by Caisse Centrale de Réassurance (CCR).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: Fault lines from CHARM database (BRGM).

Acknowledgments:
We would like to thank the industrial and academic sponsors of the RING-GOCAD Consortium managed by ASGA for their support, especially CCR. The authors also want to thanks the anonymous reviewers for their meaningful reviews that helped considerably to improve this article.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

PSHA
Probabilistic Seismic Hazard Assessment FMD Frequency-Magnitude Distribution PMD Proportion-Magnitude Distribution SIHex Sismicité Instrumental de l'Hexagone [35] FCAT-17 French CATalogue 2017 [33] Appendix A. Parameters of Rupture Plane Each region is characterized by rupture plane's parameters defined by range of values. These values are defined according to seismic and tectonic data (seismicity, faults and focal mechanisms). In this study, we use the European [49] and World [50] focal mechanism catalogues. Table A1 summaries these ranges of values for depth, azimuth, dip and movement. Since SHARE regions are large (Figure 4), range of values are wide: up to 25 km for depth, 120 • for azimuth and 40 • for dip (Table A1). All the movements (normal, strike-slip and reverse) are allowed in the stable continental region due to its wideness. Moreover, every azimuth value can be drawn in this region. Finally, since no focal mechanism is available in the oceanic crust region, all the dip (20 to 90 • ) and azimuth (0 to 359 • ) values are allowed.
Azimuth can also be described through pdf of values by deriving fault map (Figure 6a). However, these data give information only on orientation (between 0 and 180 • ) and not azimuth (between 0 and 360 • ). We thus propose to calculate pdf of orientation from fault map. Then, when an orientation O • is drawn in one of these pdf, the azimuth is chosen as equal to O • or O + 180 • . Figure A1 illustrates pdf of orientation derived from fault map for each SHARE region (except the oceanic crust one). shocks as function of magnitude. Obtained b-values for each region are summarized Table A3. Table A3. b-values computed by regressing GR law (Equation (10)) on FMD of main shocks observed in instrumental seismicity for each SHARE region. Magnitude's ranges and number of data used for regression are also detailed. Region numbers refer to Figure 4. One can see that the two extreme b-values, 0.72 and 1.36 (Table A3), are obtained with the lower number of data: 40 and 301 respectively. Thus, these values must be used carefully.