Why Are There No Earthquakes in the Intracratonic Paris Basin? Insights from Flexural Models

: Comparing nearby areas with contrasted seismicity distributions like the French Variscan Armorican Massif (AM) and the surrounding intracratonic Paris Basin (PB) can help deciphering which parameters control the occurrence or absence of diffuse, intraplate seismicity. In this paper, we examine how lithosphere temperature, fluid pressure, and frictional strength variations, combined with horizontal and bending stresses, may condition brittle, ductile or elastic behaviours of the crust in the AM and PB. We compute yield stress envelopes (YSE) and lithospheric flexure across a 1000 km-long SW–NE profile crossing the AM and PB approximately parallel to the direction of the minimum horizontal stress. Flexural models slightly better fit measured Bouguer gravity data if we apply two vertical loads on the AM and PB, with values (positive downward) ranging between − 3 and − 2.10 12 , and between 4 and 6.10 12 N.m -2 , respectively, depending on the chosen crustal composition. Our results evidence that whatever the crustal composition, bending stresses and heat flow variations alone are not sufficient to explain the difference in seismogenic behaviour between the AM and the PB. Variations in friction coefficient, in the range of standard values, are not totally satisfying either, since they do not restrain the brittle crustal thickness in the PB to less than 10 km, which is still large enough to be the locus of shallow earthquakes. Oppositely, increasing the cohesion from 10 to 80 MPa has a stronger effect on the thickness of the brittle upper crust, decreasing it from 10 to 15 km beneath the AM to 0–5 km beneath the PB. This suggests that the Mesozoic sedimentary pile can act as a sticky layer holding together basement rocks of the PB, which is equivalent to an increase in cohesion, and protects them from failure.


Introduction
Understanding the origin of seismicity in slow deforming areas requires both to decipher the stresses at the origin of earthquakes (far-field stresses due to tectonic plate motions, local stresses due to erosion, sedimentation and flexure, deep-seated processes, or any combination of these) and to characterize the mechanical response of the lithosphere to these stresses. Several hypotheses have been proposed to explain the origin of seismicity in slow deforming areas, among which are the viscoelastic behaviour of a weak lower crust responding to stress perturbations [1], stress concentration due to lithosphere heterogeneities [2], stress perturbation due to glacio-isostatic effects [3], thermal weakening of the lithosphere [4,5], or fluid flow in fault zones [6].
A classic example of enigmatic intraplate seismicity is the New Madrid Seismic Zone (NMSZ) in Central USA, where several large magnitude earthquakes (Mw > 7) have occurred, though GPS velocities depict almost no relative horizontal motions across this zone [1,5,7,8]. The Charlevoix Seismic Zone (CSZ) in Eastern Canada is another example of intraplate seismicity with almost no detectable horizontal strain [2,4,9]. With much lower magnitudes, West and Central France are also prone to diffuse intraplate seismicity, in the Variscan Armorican Massif and Central Massif (AM and CM, respectively) [10,11]. In these areas, GPS velocities also evidence null relative horizontal motions [12]. All these areas also have in common that the seismicity is at least partially localized on inherited basement faults: The Cambrian Reelfoot Rift in the NMSZ [13,14], Precambrian to Early Paleozoic faults in the CSZ [2], and Variscan strike-slip shear zones in the AM and CM [11,15], suggesting that damage zones associated with these old faults are still weaker than the surrounding crustal domains.
Understanding why there are earthquakes in a given region could also be addressed by determining why there are no earthquakes in surrounding areas, because the number of parameters (geology, heat flow, far field or local stresses) that really differ between an active seismic zone and its immediate surroundings ought to be limited. In this paper, we focus on the contrasted seismicity distribution between the French Armorican Massif and the surrounding intracratonic Paris Basin (PB). We examine how lithosphere temperature, fluid pressure, and friction coefficient variations, combined with horizontal and bending stresses, may condition brittle or elastic behaviours of the upper crust.

Geologic Setting
The AM and CM are outcropping remnants of the Paleozoic Variscan belt [16]. They are made of an assemblage of Proterozoic and Early Paleozoic metamorphic rocks of various origins: Sediments, granites, mafic, and ultramafic slivers, intruded by syn-and post-orogenic granites and overlain by late Carboniferous late-and post-orogenic basins. In the Armorican massif, the most visible Variscan structures are major Devonian to Carboniferous strike-slip crustal structures like the South Armorican Shear zone [17]. In contrast, Variscan structures of the Central Massif and in the PB are mostly flat-lying, slightly folded lower crustal thrusts, except for some late Carboniferous strikeslip faults [18].
East of the AM and North of the CM, the ~400km-wide, almost circular Paris Basin is a large intracratonic basin that developed between the Permian and the late Cenozoic [19,20]. Mesozoic sediments of the PB lie with an angular unconformity over Variscan structures of the AM and CM. To the East, these sediments are slightly tilted westwards, because of the Late Cenozoic uplift of the Vosges-Black Forest bulge [21]. The maximum thickness of Mesozoic and Paleozoic sediments is found near Paris and reaches approximately 3.5 km [20]. No major faults are visible at the surface in the PB, except the Bray-Vittel fault which is a reactivated Variscan structure ( Figure 1). Seismicity in the Armorican massif [22,23] has a diffuse character compared to interplate seismicity ( Figure 1). Some earthquake epicentres seem to cluster close to major Variscan structures, but many others cannot be attributed to any identified fault trace. Seismicity is abundant in the southern part of the Armorican Massif and on the northern border of the Central Massif. Some earthquakes occur between these two areas, at the southern extremity of the Paris basin, suggesting that active structures of these two massifs somehow connect to each other beneath the Mesozoic sediments of the basin. There is a striking difference in the epicentre density distribution between the PB and the two crystalline massifs, with much less earthquakes in the PB. Moreover, the density of earthquakes decreases rapidly with the thickness of Mesozoic-Cenozoic sediments overlying the Variscan basement ( Figure 2). The same pattern is observed further south in the Aquitaine basin [23].
Focal depths available in the seismicity catalogue [23] are probably not well constrained, owing to the low magnitude of most seismic events (≤3.0). According to [24], most hypocentres in the AM and CM are located between 0 and 10 km, i.e., in the upper crust. Earthquake focal mechanisms depict mostly normal or normal-oblique faulting [25] and the resulting stress field shows a minimum extensional stress trending ~NE-SW in the AM and PB, while the stress regime varies from extensional in the AM to strike-slip in the PB ( Figure 1) [24,26,27].  [24], [26] and [28]. Bottom: Bouguer gravity (Bureau Gravimétrique International Database), seismicity [23], and isopachs of the Paris Basin in meters (Bureau des Recherches Géologiques et Minières, Geological map of France scale 1:1,000,000). Solid line corresponds to the modelled profile.
However, according to [28], in the CM and AM the relative magnitudes of σ1 and σ2 are similar so there is no large difference, in terms of stress magnitudes, between the two stress regimes. The Bouguer gravity anomaly is typical of a stable continental lithosphere in a passive margin setting, with moderately negative values in the continental domain (−60 to −80 mGal) increasing to strongly positive ones (>150 mGal) in the oceanic domain (SW). Some short-wavelength negative and positive Bouguer anomalies in the AM and PB probably reflect density variations in the crust [11] (Figure 1).

What Can Explain the Contrasted Seismicity Distribution between PB and AM?
The main physical parameters that can increase or decrease the abundance of seismicity in a given area are the geotherm; crustal composition; orientation and friction coefficient of inherited structures; magnitude of horizontal stresses; and bending stresses due to flexure, cohesion, and pore fluid pressure in the crust. In the PB, AM, and CM, potential sources of horizontal stress are the Pyrenees, the Alps, and the Atlantic mid-oceanic ridge [29]. As there is no major actively deforming boundary between the AM, PB, and CM, we assume that the lithospheric stress does not change much between these regions. Surface heat flow is rather low in the AM and in the northern part of the PB (40-70 mW.m -2 ) and increases to more than 110 mW.m -2 southwards, towards the CM [30]. Thus, there is no striking correlation between surface heat flow and the density of earthquakes.
Except for the presence or absence of Meso-Cenozoic sediments (mainly limestones, chalk, and sandstones), there is no reason either to suppose that the crustal composition or structure differs between the PB and the AM in a way that is sufficient to explain the large difference in seismicity distribution. On the other hand, the clear drop in earthquake density with the increasing thickness of sediments in the PB suggest that the latter play a role either directly or indirectly i) by protecting underlying basement faults from alteration and subsequent development of low-friction gouge fault minerals; ii) because basement faults are at larger temperatures (60-100 °C), hence more prone to recrystallisation and healing beneath the basin than when they outcrop at the surface; or iii) because the pile of consolidated and almost undeformed Meso-Cenozoic sediments act as a sticky patch on the more intensely fractured basement, increasing its effective cohesion and preventing it from failure.

Method
We computed the yield stress envelopes (YSE) and lithospheric flexure across a 1000 km-long SW-NE profile crossing the AM and PB approximately parallel to the direction of the minimum horizontal stress, i.e., parallel to the extension direction [27]. Surface heat flow data along the profile were extracted from an interpolated grid of heat flow measurements [30]. Temperatures on the left part of the profile, in the offshore domain, were extrapolated since we do not have surface heat flow measurements.
The 2D geotherm along the profile was computed using the steady-state heat equation and assuming only vertical heat conduction ( Figure 3): with k the thermal conductivity (assumed constant), hr the characteristic depth of radiogenic heat production in the crust, and P0(x) = hrρH0(x) is the surface heat flux coming from radiogenic heat production in the crust (see Table 1 for values and units). We assume that the heat flux coming from the mantle (Qm) at the base of the crust is constant, so we can estimate the variation of the heat production from the surface heat flux Q0: (2)

Brittle Domains
The transition from elastic to brittle behaviour obeys the Mohr-Coulomb rupture criterion: Where τ is the tangential stress on a given plane, σn is the normal stress, μ is the friction coefficient and C0 is the cohesion. If we consider a medium with no weak faults, the friction coefficient is ~ 0.6 [31]. We also tested a friction coefficient of 0.3 in the crust, which simulates the presence of weaker faults. The differential brittle yielding stress σ1-σ3 is computed using the Mohr circle construction, such as with the angle between the normal to the potential slip plane and σ1. For a compressional stress regime, σ 3 is vertical and equal to the overburden pressure, whereas for an extensional stress regime σ 1 equals the overburden pressure. Pore fluid pressure favours brittle deformation by reducing the normal stress applied to the slip plane. We assume two end-member cases with either a dry crust (i.e., with no fluid pressure) or a "wet" crust with hydrostatic fluid pressure. We do not consider overpressured fluids the AM as there is no cap-rock lying upon basement aquifers that could help building-up fluid overpressure. We also consider a wet crust (i.e., hydrostatic fluid pressure) in the PB as, even if most aquifers are contained within the sedimentary layers, it is unlikely that the basement is completely dry.

Ductile Domains
The transition from elastic to ductile behaviour is computed using a temperature and stressdependent dislocation creep power law: where ε is the strain rate (assumed constant and equal to 1e -17 s -1 ), n is a power law exponent, H is the creep activation energy, R is the gas constant, σ1-σ3 is the differential yielding stress, and T is the temperature in Kelvins (Table 1). We used a dry olivine rheology for the upper mantle, and either diorite or diabase rheologies for the crust [32,33]. The YSE is the contour of the minimum yielding stress, i.e., the differential stress which conditions the transition between elastic and either brittle or ductile behaviours at each depth. It is then used to compute the effective elastic thickness (EET) of the lithosphere.

Effective Elastic Thickness and Flexure
We consider as pertaining to elastic crust and mantle lids the domains where the yielding stress (either in ductile or brittle domains) is larger than a threshold of 20 MPa [34]. If the crust and mantle are mechanically coupled (i.e., there is no weak ductile lower crust), EET is simply the sum of the crust and mantle elastic layer thicknesses, hcrust and hmantle, respectively. If there is a mechanical decoupling between the upper crust and upper mantle, then The effective elastic thickness is then used to compute the flexural rigidity D and solve the 1D flexure equation using a 1D finite element model: where E and ν are Young's modulus and Poisson ratio, ρm and ρf are the mantle and infill (water offshore and air onshore) densities, w is the plate deflection, and q is the vertical load applied to the plate. The latter comprises both the present-day positive topography load and possible additional forces due either to dynamic processes in the mantle, erosion at the surface, or to buried loads in the lithosphere. Plate deflection cause bending stresses σb that modify the thickness of the elastic layers, hence the effective elastic thickness. In addition, horizontal stresses σh applied to the lithosphere also alter the YSE and modify EET [34]. The resulting differential stress resolved on the lithosphere is σd = min(|YSE|, |σh+σb|). We iteratively compute EET and w until the model converges towards a stable solution, while imposing null deflection at both plate extremities. The Moho deflection is then used to compute a synthetic Bouguer gravity anomaly (assuming a crustal density ρc), which is then corrected for the sediment infill of the PB and compared to Bouguer gravity data.

Results
Whatever the chosen crustal composition, we find that flexural models slightly better fit Bouguer gravity data if we apply two vertical loads (positive downward), at x = 400 km and x = 750 km, with magnitudes ranging between −3 and −2.10 12 , and between 4 and 6.10 12 N·m -2 , respectively, depending on the chosen crustal composition (Figures 4 and 5). The first load could correspond to isostatic uplift of the AM in response to erosion or post-glacial adjustment, which has already been suggested as a possible source of stresses in the AM [3]. The second one might correspond to a buried load beneath the deepest part of the sedimentary basin, which may explain why this basin experienced subsidence. A positive value of 5.10 12 N.m -2 corresponds for example to a positive density perturbation of 100 kg.m -3 on a 1 km × 1 km × 5 km crustal block. In the following models, we keep these loads in order to insure a good fit to gravity data (see Discussion). Moreover, as bending stresses due to flexure modify the yield stress envelopes, it is important to test whether these additional loads, although small, could explain the seismogenic behaviour of the AM and PB.  Table 1 for details). Flexural models with a strong continental crust (dry diabase), a friction coefficient of 0.6, a cohesion of 10 MPa, and no horizontal stresses predict that the continental crust behaves elastically along the largest part of the profile, except in places where bending stresses are large enough to reach the yielding stress of the ductile lower crust. The elastic cores of the crust and upper mantle are coupled except around the points where bending stresses are the highest, which causes drastic drops in EET from > 80 km (coupled case) to 25 km (decoupled case) along the profile ( Figure 6). In this first model, differential stresses are much too low to cause any permanent deformation in the lithosphere.
If a horizontal extensional stress of 150 MPa is added, all the uppermost crust falls in the brittle field and there is no striking difference between the rheological layering of the AM and PB ( Figure  7). Compared to the previous test, in this case the elastic crustal and mechanical layers are always mechanically decoupled, so that EET varies between 20 km and 35 km.  If we use a weaker crust rheology (diorite) with a hydrostatic pore fluid pressure and a reduced friction coefficient of 0.3, we obtain results that are not satisfying: Without horizontal extension, there is only a very thin (0-3 km) and discontinuous brittle crustal layer, then a ~10 km-thick elastic lid, below which the crust behaves as a weak ~20 km-thick ductile layer. The EET varies between 20 and 30 km, with some higher peaks where plate curvature is null (Figure 8). The extremely thin (or inexistent) brittle crust and the thick ductile lower crust predicted by this model are difficult to reconcile with the abundance of earthquakes in the AM. Moreover, there is still no visible difference in the mechanical behaviour of the crust between the AM and the PB.  For the same crustal rheology and with a low horizontal extensional stress of 50 MPa, a thicker and more continuous part of the uppermost crust (~10 km) falls in the brittle field, the crustal elastic lid is reduced to less than 2 km, and the rest of the crust is weak and ductile. Moreover, here again there is no difference of crustal mechanical behaviour between the AM and the PB (Figure 9). This first set of models evidence that if we keep a constant crustal rheology along the profile, neither bending stresses nor heat flow variations can explain the contrasted seismogenic behaviour between the AM and PB. We then consider two different crustal domains. Multiple tests allow us to determine that, overall, a "strong" crust composition (diabase) better suits than a weak (diorite) one for the whole profile. In a next series of models, we test the possibility of a different friction coefficient and/or a different cohesion between the left part of the profile (from 0 to 540 km), which corresponds to the Armorican Massif domain, and the right part of the profile, corresponding to the Paris Basin. We keep a constant horizontal extensional stress of 100 MPa and a hydrostatic fluid pressure for the whole profile. Varying the friction coefficient allows us to simulate the effect of different material properties in the crust (due to variable alteration or recrystallization of inherited fractures for instance), whereas the cohesion is modified in order to simulate the presence or absence of a sticky patch of consolidated sediments overlying the basement.   In the model with two different friction coefficients of 0.3 and 0.7 for the AM and PB, respectively, the thickness of the brittle crust decreases from ~ 20-25 to 10-12 km between the left and right parts of the profile (Figure 10). In the right part of the profile (corresponding to the PB) the elastic crust is much thicker than on the left part (>20 km) but the brittle crust is thick enough to be the locus of numerous shallow earthquakes, which is not consistent with observations. As the crust and mantle are always decoupled, the EET varies between 20 and 40 km. Finally, if we assume a change in cohesion (10 to 80 MPa) between the left and right parts of the profile without changing the friction coefficient (0.6), we can adequately reduce the brittle crustal thickness to 0-5 km and increase the thickness of the crustal elastic lid in the PB (Figure 11), while the AM crust is characterized by three approximately equal brittle, elastic and ductile layers. In this model, we obtain a large contrast in the thickness of the brittle crust between the AM and PB, which seems more consistent with earthquake distribution. The effective elastic thickness varies between 20 and 40 km and the crust and uppermost mantle are decoupled.

Discussion and Conclusion
Our results evidence that whatever the chosen crustal composition, bending stresses and heat flow variations alone are not sufficient to explain the difference in seismogenic behaviour between the AM and the PB. Variations in friction coefficient between AM and PB, in the range of tested values (0.3 to 0.7, respectively) are not totally satisfying since they do not restrain the brittle crustal thickness in the PB to less than 10 km, which is still large enough to be the locus of shallow earthquakes. As a friction coefficient of 0.7 is in the upper range of values commonly inferred for intact crustal rocks, we doubt that a much larger (i.e., >0.8) value could realistically be used. In fact, a friction coefficient of 0.9 still increases the thickness of the elastic lid, but it does not sufficiently reduce the thickness of the brittle crust in the PB to explain the lack of earthquakes. Oppositely, increasing the cohesion from 10 to 80 MPa has a stronger effect of the thickness of the brittle upper crust, decreasing it from 10 to 15 km beneath the AM to less 0-5 km beneath the PB. A cohesion of 80 MPa is in the upper bound of expected values for strong, cohesive rocks [35]. Our results therefore suggest that the unfaulted Meso-Cenozoic sedimentary pile acts as a sticky layer holding together the faulted basement rocks of the PB. This is equivalent to an increase in cohesion ( Figure 12). Mechanical coupling in the PB between basement rocks and the overlying sedimentary pile has already been proposed [36] and strengthens this interpretation.
Another outcome of our models is the use of negative and positive vertical loads on the central parts of the AM and PB, respectively. Since the fit to Bouguer gravity data is only slightly better with these line loads than without ( Figure 5), we cannot claim that they are absolutely necessary. Flexural models with the same setting as Figure 11 but without line loads produce a large contrast in brittle crustal thickness between the AM and PB. Still, it was important to estimate the maximum magnitude of these vertical loads and our study suggests that they cannot be responsible for the contrasted seismicity distribution between the AM and PB. However, they can explain, for instance, i) the longterm subsidence of the PB; ii) the relative uplift of the AM [37]; and iii) the supposed activity of the reverse Bray fault north of the PB [38].
Our thermal model with a constant heat conductivity could also be questioned, in particular because of the presence of low conductivity Triassic evaporites in the PB [20]. Indeed, temperatures computed with a thin (1 km) low-conductivity (0.5 W·m -1 ·K -1 ) layer are up to 100 K larger than in our model. As a result, ductile-elastic transitions in the crust and mantle beneath the PB are uplifted by 7 and 9 km, respectively. However, this thermal blanketing has no effect on the depth of brittle-elastic transitions, which is not temperature dependent, hence on the brittle behaviour of the crust.
Unlike in previously published studies [39], [40], we cannot use here the depth distribution of earthquakes to better constrain the crust rheology. Due to the low magnitude of events and relative looseness of the seismic network, we can only assume that most earthquakes are contained within the first ~10-15 km of the crust.
The differences in the seismic behaviours between the AM and the PB come from the pile of the sedimentary layers on top of the PB. We here show that these sediments may act as a sticky patch, which is equivalent to an increase in cohesion. This strength increase prevents the occurrence of seismicity as long as there are no stress perturbations. However, it does not mean that the fault state is far from failure. Indeed, transient perturbations of local stress, fault strength, or fluid pressure may generate earthquakes in those area [41], as seen in, e.g., the Oklahoma state where a small fluid pressure increase lead to large earthquakes and cascading seismicity [42].
Funding: This research received no external funding.