# Exceptional ground motion during the shallow Mw 4.9 2019 Le Teil earthquake, France

Mathieu Causse, Cécile Cornou, Emeline Maufroy, Jean Robert Grasso, Laurent Baillet, Elias El Haber

► To cite this version:

Mathieu Causse, Cécile Cornou, Emeline Maufroy, Jean Robert Grasso, Laurent Baillet, et al.. Exceptional ground motion during the shallow Mw 4.9 2019 Le Teil earthquake, France. Communications Earth & Environment, 2021, 2 (14), 10.1038/s43247-020-00089-0. http://hal-03123311

HAL Id: hal-03123311

https://hal.science/hal-03123311

Submitted on 31 Jan 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for updates

ARTICLE

https://doi.org/10.1038/s43247-020-00089-0

OPEN

# Exceptional ground motion during the shallow Mw 4.9 2019 Le Teil earthquake, France

Mathieu Causse<sup>D 1✉</sup>, Cécile Cornou<sup>1</sup>, Emeline Maufroy<sup>D 1</sup>, Jean-Robert Grasso<sup>D 1</sup>, Laurent Baillet<sup>1</sup> & Elias El Haber<sup>1</sup>

An unusually damaging Mw 4.9 earthquake occurred on November 11, 2019 in the south east of France within the lower Rhône river valley, an industrial region that hosts several operating nuclear power plants. The hypocentre of this event occurred at an exceptionally shallow depth of about 1 km. Here we use far-field seismological observations to demonstrate that the rupture properties are consistent with those commonly observed for large deeper earthquakes. In the absence of strong motion sensors in the fault vicinity, we perform numerical predictions of the ground acceleration on a virtual array of near-fault stations. These predictions are in agreement with independent quantitative estimations of ground acceleration from in-situ observations of displaced objects. Both numerical and in-situ analyses converge toward estimates of an exceptional level of ground acceleration in the fault vicinity, that locally exceeded gravity, and explain the unexpectedly significant damage.

<footnote n="1">¹ISTerre, Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, Univ. Gustave Eiffel, 38000 Grenoble, France. ✉email: mathieu.causse@univ-grenoble-alpes.fr</footnote>
Earthquake ground motion results from a sudden release of elastic energy, accumulated over decades to millennia, in the shape of a slip distribution along a fault. Surface ground motion transmits forces to structures and thus, predicting the ground motion is the basis of earthquake risk assessment. Strong motion records are typically used by the earthquake engineering community to design key structures like power plants, bridges or dams, and to analyze the response of existing buildings for a potential future earthquake. One major issue is that near-fault strong motion recordings are still too scarce to fully understand the physics involved in the highly variable ground motion within a few kilometers from the rupture plane (e.g., ref. ¹,²). Consequently, ground motion prediction models are not calibrated for near-field motion and are accurate only at fault-plane distances greater than ~10 km (ref. ³). Numerical ground motion simulation techniques, including physical description of the fault rupture process, have then become popular to investigate near-source characteristics of ground motion (e.g., ref. ⁴).

On November 11, 2019, a M<sub>w</sub> 4.9 earthquake occurred in the south east of France (referred to as the Le Teil event, Fig. 1a) close to the city of Montélimar within the lower Rhône river valley, a moderate seismicity area. This industrial region hosts several operating nuclear power plants. The unexpected level of damage for this magnitude (epicentral intensity VIII, ref. ⁵) attests an intense ground shaking. Unfortunately strong motion was not recorded by accelerometers near the fault, in the damage area. The earthquake causative fault, which was not considered as active, is part of the Cévennes fault system (ref. ⁶,⁷). Seismological, geodetic and field observations indicate a reverse-fault mechanism, a rupture area of about 4 km by 1.5 km, a hypocenter at about 1 km and surface rupture evidences (ref. ⁷,⁸). Such a shallow location is exceptional. Moderate-sized earthquakes overwhelmingly occur at seismogenic depths larger than ~5 km in stable continental regions (ref. ⁹). The shallow location may be related to quarry extraction directly above the rupture area (ref. ¹⁰⁻¹²). In a few cases, small to moderate-sized anthropogenic earthquakes with reverse-faulting style are reported in the immediate proximity of quarries (e.g., ref. ¹³⁻¹⁶). Superficial hypocenter depth is reported as a common feature for anthropogenic seismicity (e.g., ref. ¹⁷,¹⁸).

Near-field seismological observations of moderate-sized superficial events (depth <3 km) are inexistent in France and only a very few are reported worldwide. The European strong motion database (ref. ³) contains only eight recordings (over 23,014) of superficial events recorded <5 km from the rupture plane, all of them being from earthquakes with local magnitude lower than 4.5. Two cases of superficial events with M<sub>w</sub>~5 and causing exceptional ground motions (i.e., with Peak Ground Acceleration—PGA—exceeding 5 m/s²) were reported in California (ref. ¹⁹) and in Nevada (ref. ²⁰).

In the absence of near-fault strong motion recording for the Le Teil earthquake we conduct numerical predictions of the ground acceleration (Fig. 1b). First we perform a source analysis and generate a suite of rupture models, describing the space-time distribution of slip across the fault. The rupture models are calibrated based on the Interferometric SAR slip model proposed by (ref. ⁸). Next we characterize the 1D structure of the earth crust using seismic noise recorded at temporary seismological stations installed after the earthquake in the fault vicinity (ref. ⁸). Finally we obtain a set of ground acceleration time histories on a dense array of virtual stations, centered on the fault. Furthermore, we collect in-situ observations providing independent estimates of the ground acceleration, such as displacements of funeral slabs in cemeteries and displaced rocks above the fault. Both numerical analysis and in-situ measurements of displaced objects converge toward exceptional ground acceleration, exceeding gravity in the immediate fault vicinity. Our study provides new insights on seismic hazard due to superficial seismicity in stable continental regions.

## Results

**Source analysis and rupture modeling.** Defining proper rupture models is a crucial point because near-fault surface ground motion is sensitive to the details of the fault slip history. Theoretical and numerical studies indeed show that near-fault ground acceleration is primarily controlled by local processes on the fault, occurring near the recording site, including local stress drop, strong spatial variations of the rupture velocity (in particular

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

**Fig. 1 Methodology to compute ground acceleration.** a Map of the Le Teil area. The colored rectangle shows one realization of rupture model (chosen within the 2000 generated models), with final slip indicated by the color scale (max = 40 cm in red), and white contour lines indicating the rupture front every 0.5 s. The rupture plane reaches the surface and dips 58° toward the southeast. The star denotes the hypocenter. Modeling of funeral slab motion is performed at cemetery C1, including a large majority of the observations (Table S3, Fig. S7). We observed only a few (<5) displaced slabs at C2 and C3, while all slabs remained stable at C4 and C5. The orange polygon is the area where displaced rocks are observed. b Flow chart describing the strategy to compute near-fault ground acceleration.
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Velocity model. (Left panel) Shear-wave velocity profile (black line indicates the best misfit profile, gray lines indicate the ensemble of inverted $V_s$ profiles that explain dispersion data within their uncertainty bound) overlaying the geological cross-section modified from (ref. ⁷) (licensed under CC BY 4.0.). (Right panel) Pictures of the outcropping Hauterivian and Valanginian marls and limetones.</figcaption></figure>

when the rupture front hits the fault edges) and local fault slip acceleration (e.g., refs. ²¹⁻²³). We use a kinematic representation of the rupture process. Six properties must then be defined, namely the rupture size, the final slip distribution, the rupture duration, the position of the rupture initiation, the rupture velocity and the local slip velocity. As shown hereafter, most of them are constrained by observations, except the slip velocity that needs to be calibrated based on past earthquake studies.

Thanks to the shallow rupture, the co-seismic displacement field could be mapped by Interferometric SAR and inverted to provide a final slip distribution (ref. ⁸), that we use to constrain the low frequency slip features (Methods). Slip shows a reverse-fault mechanism and is characterized by two main slip patches localized between 0 and ~1 km depth on a roughly 4 km*1.5 km fault plane dipping 58°SE, resulting in a magnitude $M_w$ = 4.9. Analysis of the farfield seismological recordings gives the same magnitude (moment tensor inversion, ref. ⁷), indicating that the slip obtained from InSAR is primarily seismic.

In order to constrain the rupture duration, we use seismograms recorded by the real-time broadband and accelerometric permanent stations of the French seismologic and geodetic network (RESIF) (59 stations with source-site distances from 27 to 290 km, Fig. S1a). We employ a spectral approach and perform a Bayesian estimation of the corner frequency (Methods, Fig. S1b). The later is related to the rupture duration given a source model. Our analysis yields rupture duration of 2.0 s (confidence interval at 68%, 68% CI: 1.4–3.4 s). The Brune (1970) stress drop²⁴ we obtain is 1.0 MPa (68% CI: 0.2–3.4 MPa). This value is close to the average of worldwide reported stress drops (e.g., refs. ²⁵⁻²⁷). Furthermore, the azimuthal distribution of the apparent rupture duration (i.e., “seen” from any station) is classically used to track the main rupture direction (e.g., refs. ²⁸⁻³⁰). The corner frequency shows weak azimuthal dependence, revealing a bilateral rupture, i.e., initiation about midway between the SW and NE rupture terminations (Methods, Fig. S2). This is consistent with the analysis by (ref. 11) and with the hypocenter location available at the time, based on inversion of farfield seismological data and placing the hypocenter at the fault center at 1 km depth with spatial uncertainty of 1 km (ref. ⁷). The rupture may have started slightly southwestward as suggested by coherency analysis of high-frequency waveforms at two neighboring farfield stations

(ref. ³¹). Knowing the position of the rupture initiation, the rupture duration and the rupture dimension, we infer an average rupture velocity of ~1.8 km/s. Considering the 1D velocity model developed in this study (section Velocity model), it is between 50% and 90% of the shear-wave velocity ($V_s$), as reported for most large earthquakes (e.g., refs. ³²,³³).

Furthermore, it is necessary to define how local slip evolves with time (slip velocity function). We use a function in which slip starts with a short acceleration phase and then slowly decelerates until the final slip is reached (ref. ³⁴), as shown by laboratory experiments and dynamic rupture simulations (e.g., refs. ³⁴,³⁵). We impose constant slip duration of 0.3 s (68% CI: 0.2–0.4 s) as proposed by³⁶ for a $M_w$4.9 event. The average slip velocity over the fault is then 0.5 m/s, as commonly reported from near-fault seismological data (e.g., refs. ³²,³⁷) and observations during fault scarp formation (e.g., refs. ³⁸,³⁹). The peak slip velocity locally reaches 4 m/s (Fig. S3a, ref. ⁴⁰). The duration of the slip acceleration phase is poorly constrained, but we checked that the ground acceleration time series presented hereafter are not sensitive to this parameter in the considered frequency range (<5 Hz).

In order to account for source heterogeneity involved in the high-frequency seismic energy (e.g., ref. ⁴¹), we introduce small-scale variability of the rupture parameters by considering self similar spatial distributions. Self-similarity of slip distribution is supported by observations of slip images obtained from geodetic or seismological data and topography analysis of exhumed faults (ref. ⁴²). Finally, we obtain 2000 rupture models representing the uncertainties in the above-defined rupture parameters, described by probability density functions (Methods, Fig. S3).

**Velocity model.** We derive a 1D velocity model of the earth crust based on three-component beamforming method (ref. ⁴³) using the seismic noise recorded at temporary post-seismic stations installed in the fault vicinity (Fig. 1a) and the Conditional Neighborhood Inversion Algorithm (ref. ⁴⁴) (Methods, Fig. S5). The shear-wave ground velocity profile (Fig. 2) exhibits materials with increasing stiffness from the surface to 1.2 km depth over-laying less competent deposits. This peculiar ground velocity profile with the presence at depth of softer material is consistent with the geological settings of the area (ref. ⁴⁵) and deep boreholes
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Numerical prediction of ground acceleration. Top: simulated median values of the Peak Ground Acceleration (PGA) between 0 and 5 Hz on the array of virtual stations. The inter-stations distance is 500 m. Bottom: PGA variability, computed over 2000 realizations of rupture models, expressed in terms of the standard deviation of the log10 values. The dashed rectangle represents the surface rupture projection. The thick black polygon line is the area where displaced rocks are observed.</figcaption></figure>

in the region (geological logs BSS002ARWX, BSS002ASXR, BSS002ASEZ available at http://infoterre.brgm.fr) that show, from surface to depth, alternating layers of marls and limestone from the lower Cretaceous, a more competent thick Jurassic limestone layer, and thick claystones from lower to upper Jurassic age that overlay a Triassic bedrock. At the vicinity of the fault zone, the geological cross-section proposed by (ref. ⁷) reports a 0.7 km thick Barremian limestone layer overlaying a 1.7 km thick formation composed of marls and limestones from Hauterivian and Valanginian age, resting on Tithonian formation. Outcropping Hauterivian and Valanginian formations (Fig. 2) indicate that the Valanginian formation is particularly rich in soft marls compared to Hauterivian formation, consistently with observations in the deep regional boreholes and in the Vonvocian basin (ref. ⁴⁶). The 1D velocity model is thus consistent with geological information leading to the following interpretation: 0.6 km thick Barremian limestones with Vs ~2 km/s and 0.6 km thick Hauterivian marls and limestones with Vs ~3.6 km/s overlaying Valanginian marls with Vs ranging from ~1.2 km/s to 2.3 km/s. The velocity contrast at about 2.2 km depth would correspond to the Tithonian formation. The rupture is confined between the surface and the soft Valanginian marl layer, which probably acted as a barrier to the rupture propagation. It likely initiated in the stiff Hauterivian formations, expected to be more brittle, which is consistent with the hypocenter position at 1 km depth inferred from farfield seismological data (ref. ⁷).

**Numerical prediction of ground acceleration.** In order to simulate ground acceleration, the seismic energy radiated by each rupture model is propagated in the inferred best misfit velocity model toward the array of near-fault virtual stations. The array is composed of 234 stations with interstation distance of 500 m and is centered on the fault (Fig. 3). Then we compute the impulse

response of the medium and following the Representation theorem (ref. ⁴⁷), we convolve the medium response with the slip velocity and obtain 2000 ground motion time series up to 5 Hz for each station (Methods).

The results indicate intense ground shaking above the fault rupture projection, primarily vertical, with horizontal PGA up to 7.5 m/s² (68% CI: 3.0–19 m/s²) and vertical PGA exceeding gravity (13 m/s², 68% CI: 5.8–29 m/s²) above the main slip area (Figs. 1a and 3). The decay of ground acceleration with distance to rupture is however very fast. In addition, directivity effects are known to cause ground motion amplification at a given site as the rupture propagates toward the site and when the direction of slip points toward the site (e.g., refs. ⁴⁸,⁴⁹). For a reverse-fault mechanism with a dip angle of ~60°, such effects are assumed to be very weak along the fault strike direction—as observed in our simulations—and maximized along the up-dip direction. Nevertheless, the weak amount of slip above the nucleation area for the Le Teil event (Fig. S3a) results in the absence of up-dip directivity effects and consequently, ground accelerations are lower toward NW.

It is known that PGA for small earthquakes can be strong without causing damage, because it may be controlled by energetic frequencies that are much higher than the typical resonance frequencies of civil structures (>20 Hz). However, our simulations result also in large values of the Peak Ground Velocity (PGV) for frequencies below 5 Hz (exceeding locally 1 m/s, Fig. S4). Moreover, the ratio PGA/PGV, which equals about 2π multiplied by the dominant frequency, has values ranging between 10 and 15 s⁻¹, implying a dominant frequency around 2 Hz and hence accelerations prone to damage structures (ref. ³), in agreement with the observed Intensity of VIII in the fault vicinity.

**Near-field observations of displaced objects.** In order to obtain independent estimations of the ground acceleration, we collected
<figure><img src="image_7.png" /><figcaption>**Fig. 4 Displaced rocks.** Examples of observed displaced rocks as upthrown stones, freshly fractured and broken at the impact on ground, imply a vertical acceleration exceeding gravity during the Le Teil event. a-d Broken stones, apparently flung laterally. For c and d, the bottom figures show the reconstruction from fragments. The detailed observations and site locations for all displaced rocks are available in Fig. S6 and Table S2. All the observed site locations, which cluster close to the SW branch of the surface rupture, are reported as a polygon on Figs. 1 and 3, for comparison with the simulated seismic slip and ground acceleration.</figcaption></figure>

in-situ observations of displaced rigid objects, either natural or man-made, that fulfill the following criteria: (1) there is evidence that displacement is due to the earthquake; (2) the objects loosely lie on the ground; (3) their structure is simple enough so that their displacement can be related to ground acceleration using simple mechanical concepts. We compiled a database of displaced rocks and a database of funeral slabs displacements.

After the Le Teil event, we observed, in specific places, that tens of stones and rocks, which are numerous in this back country of smooth hills and scattered limestone outcrops, had been displaced and sometimes freshly fractured and broken (Fig. 4, Methods, Table S2, Fig. S6). Freshness of cuts and fractures and the consistency of several tens of observations indicate that some of these rocks were tossed into the air during the earthquake, with breakage occurring at the time of impact, implying that vertical acceleration exceeded gravity (e.g., refs. 50,51). All the evidences for upthrown and displaced rocks are located in the vicinity of the SW part of the surface fault trace (Fig. 1a and S6), an area that remarkably coincides with the area where simulated vertical PGA exceeds gravity (Fig. 3).

In addition, tens of funeral slabs were displaced in a cemetery (C1) located about 200 m from the northern rupture termination (Fig. 5a, Methods, Table S3, Fig. S7). The cemetery is located within an area of Intensity VIII on the EMS-98 scale (ref. 52). This translates into a horizontal PGA of 2.8 m/s² (68% CI: 1.5–6 m/s²) from empirical relationships (ref. 53), which coincides with the median horizontal PGA value we simulated at this place (2.8 m/s², 68% CI: 0.5–7.9 m/s²). In order to obtain independent estimations of the ground acceleration, we model the slab motion by conducting numerical simulations based on a Coulomb friction model. The method uses the suite of 2000 three-component ground motion time series modeled in cemetery C1 as input and provides the slab displacement time series in the

EW-NS plane (Methods, Fig. 5b). The slabs can move with respect to ground if:

$$ \frac{\sqrt{a_{EW}^2(t) + a_{NS}^2(t)}}{\mu} - a_Z(t) > g \quad (1) $$

where $a_{EW}$, $a_{NS}$ and $a_Z$ represent ground acceleration in the NS, EW and vertical directions, respectively, $\mu$ is the static friction coefficient and $g$ is gravity. 65% of the simulated ground motion time series verify Eq. 1, and hence can induce a slab displacement. The slab motion is illustrated in Fig. 5b for one realization of ground acceleration that well represents the field observations. Equation 1 is verified soon after the S wave arrival, resulting in a ground displacement toward the south. This S-wave pulse induces a force on the slab in the opposite direction (north), but the headstone blocks the slab motion. From the simulations, the slab starts sliding with the arrival of a stopping phase (pulse of high-frequency ground motion) generated by the rupture arrest at the northern fault termination. This stopping phase is associated with an abrupt change of ground displacement polarity and a force rotating toward the south, coinciding with the final southward slab displacement. A global analysis (Fig. 5c) shows that a horizontal PGA of 3 m/s² generates slab motion of 0.8 cm (68% CI: 0.3–1.2 cm), while we observed in cemetery C1 a median displacement of 4.5 cm (68% CI: 2–7 cm). Such a displacement corresponds to a simulated horizontal PGA of 4.8 m/s² (68% CI: 4.3–5.4 m/s²). This indicates that the numerical ground motion simulations probably provide a lower bound of the ground acceleration that actually happened at the cemetery C1. Indeed, local site amplification effects are not considered in this study. Finally the median simulated PGA values at cemeteries C2–C5 are in the range 1.5–2 m/s², consistent with the quasi absence of observed displaced slab.
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Displaced funeral slabs. a Illustration of the typical funeral slab layout in cemetery C1 (see Figs. 1 and 3 for C1 location). The slabs back onto N70°-oriented headstones. Most of the slabs lying on the SE side were displaced during the ground shaking, while slabs lying on the other side experienced smaller displacements or did not move (Table S3, Fig. S7). Note that some of the slabs collapsed, as shown on the bottom-right inset figure. b Example of ground motion acceleration time history (top) and resulting modeled slab displacement using a friction coefficient of 0.2 (bottom) (Methods). The color scale on the bottom figure points to the different time steps as defined on the ground motion time histories, for a better visualization of the time evolution of the slab displacement. The color scale defines the time window during which Eq. 1 is verified, that is when the inertial force is sufficient to induce sliding of the slab. The arrows on the bottom figure indicate the direction and amplitude of forces transmitted to the slab at different time steps. The ground displacement (light gray) and slab displacement (dark gray) diverge at the arrival of the stopping phase generated by the rupture arrest, about 200 meters from the cemetery. The relative slab displacement is shown in purple, indicating a final displacement of 6 cm in the direction N168°. c 2D-probability density function of simulated PGA (geometrical mean of horizontal components) and modeled relative final slab displacement, obtained from 2000 realizations of rupture models. The median final relative slab displacement measured on field (4.5 cm) yields a median PGA value of 4.8 m/s² (CI 68% 4.3–5.4 m/s²) in the cemetery during the Le Teil earthquake.</figcaption></figure>

## Conclusion

The 2019 M<sub>w</sub>4.9 Le Teil earthquake occurred in a stable continental region at unusually shallow depth (~1 km). Such superficial events are rare and thus their rupture characteristics are very poorly documented. The level of ground motion that they can produce is basically unknown. We demonstrated that the average rupture properties (stress drop, rupture- and slip- velocity) are consistent with the ones commonly observed for large deeper earthquakes. We also showed that the unusually shallow rupture with generation of energetic seismic waves up to the ground surface resulted in an exceptional and destructive level of ground acceleration exceeding gravity on some localized patches in the immediate fault vicinity, as evidenced by the consistency between the numerical predictions and the body of in-situ observations. Then, our simulation results point on a fast decay of motion amplitude with fault distance for the first few kilometers. Inside the near-source area, the simulated ground motion shows a high level of heterogeneity that is fully related to the rupture characteristics along the fault plane (distribution of slip and rupture velocity, position of hypocenter).

Our results dramatically change the perception of the impacts of superficial moderate earthquakes on local hazard assessment. The scarcity of observations of Le Teil-event type implies that they are not represented in regulatory assessments of seismic hazard based on empirical ground motion models (e.g., ref. 54). Earthquake recurrence models and ground motion models that would be accurate for such superficial, possibly triggered events are not available. This raises the question of how prepared the society can be to deal with this low probability high impact- type of event. This is particularly critical in areas where long term safety of nuclear plants is involved. The Le Teil earthquake occurred in the vicinity of active nuclear power plants located within the Rhône river valley, 12 km north- and 20 km southward, respectively. The earthquake causative fault, which was not considered as active, is part of the Cévennes fault system, this later extending eastward across the Rhône valley a couple of kilometers from the northern nuclear power plant (ref. 6,7). Our analysis points on the key role of the 1.2 km thick superficial stiff rock matrix on the unusually strong ground shaking we quantify for such an earthquake size. Such geological units being common
for this region (ref. ⁴⁵), future geological and geophysical analyses should characterize the extent and thickness of the geological units and the seismogenic potential of the local faults, so as to constrain the probability of a superficial earthquake (either tectonic or anthropogenic) in the immediate vicinity of the power plants.

## Methods

**Spectral analysis of seismological data.** For each station of the RESIF network used in our analysis (Fig. S1a), we compute the seismic moment $M_0$ and the corner frequency $f_c$. Due to the unusually weak number of aftershocks, we could unfortunately not find adequate small earthquake recordings to represent the details of the wave propagation in the earth crust. We then use a simplified spectral approach (e.g., ref. ⁵⁵). We select the P wave (time window of 12 s starting from the P arrival, ref. ⁵⁶) and compute the RMS of the 3-components Fourier displacement amplitude spectra. According to²⁴, the level of the low frequency plateau $\Omega_0$ is:

$$ \Omega_0 = \frac{2M_0 U_{\phi\theta}}{4\pi\rho c^3 R} \qquad (2) $$

where $U_{\phi\theta}$ is the average radiation pattern (0.52 for P wave), c is the P wave velocity, $\rho$ is density and R is the hypocentral distance. For each station, we then normalized the displacement spectrum by $2U_{\phi\theta}/(4\pi c^3 R)$ so that the low frequency plateau provides the seismic moment. Finally, we stack all the spectra and model the resulting average spectrum as (refs. 56,57):

$$ \Omega(f) = \frac{e^{-n\eta f/Q}}{\left[1 + \left(\frac{f}{f_c}\right)^{\eta n}\right]^{1/\eta}} \qquad (3) $$

where t is the travel time of the considered wave and Q is the quality factor of the considered wave.

We then perform a Bayesian estimation of the model parameters $f_c$, $M_0$ and Q (Supplementary Methods), yielding median values $M_w = 4.85$ (68% CI: 4.65–5.05) and $f_c = 0.44$ Hz (68% CI: 0.26–0.66 Hz) (Fig. S1b). The joint posterior probability density functions $P_{f_c,m_0}(f_c, m_0)$ and $P_{Q,f_c}(Q, f_c)$ indicate that the $f_c$ uncertainty mainly arises from a strong tradeoff between Q and $f_c$ (Fig. S1c). Following⁵⁸, we obtain that the rupture duration for P-waves is between 0.84 $f_c$ and $f_c$ assuming a symmetric circular crack rupture model and a rupture velocity $V_R$ between 50% and 90% of the S-wave velocity $V_S$ (consistent with the inferred average $V_R$ of ~1.8 km/s and our velocity model), i.e., an average duration of 2.0 s (68% CI: 1.4–3.4 s).

We next compute the corner frequency for various classes of source-station azimuths, and compare the azimuthal variations with the theoretical predictions for a simple horizontal line source (Fig. S2). For a unilateral rupture:

$$ f_{cU}(\theta) = \frac{1}{T(1 - \alpha \sin(i) \cos(\theta))} \qquad (4) $$

where T is the rupture duration, $\alpha$ is the ratio $V_R/V_S$, i is the take off angle and $\theta$ is the source-station azimuth. For a bilateral rupture,

$$ f_{cB}(\theta) = \min[f_{cB}(\theta), f_{cB}(-\theta)] \qquad (5) $$

Despite the uncertainty on $f_c(\theta)$ and considering various values of i and $V_R$, the results indicate a bilateral rupture (Fig. S2).

**Kinematic rupture modeling.** We generate a suite of 2000 rupture models with the code developed by⁵⁹. The rupture area is represented by a rectangle with fixed length L and width W of 4.2 km and 1.4 km, respectively. The distributions of slip and rupture time are self similar ($k^{-2}$ decay in the wavenumber domain), and the local slip duration (rise time) is constant over the fault plane. The rupture duration follows a lognormal distribution with median of 2 s and standard deviation of 0.15 (in log10 units) (Fig. S3d). The along-strike position of the rupture initiation follows a normal distribution with median 0.5*L and standard deviation 0.1*L, while the along-dip position is uniformly distributed between 0.5 and 1 km depth. The values of the local slip duration are uniformly distributed between 0.15 and 0.45 s. The method ensures that source functions are compatible with a “ω$-2$” model²⁴ (Fig. S3d). For a specific application to Le Teil earthquake, we bring the following modifications to the code: (1) the self similar slip distributions, which have random phases, are low pass filtered below a wavenumber of $1.5 e^{-3} m^{-1}$ and replaced with the InSAR slip, high pass filtered above the same wavenumber; (2) the slip velocity function is defined according to (ref. ³⁴); (3) the rupture time perturbations are replaced by rupture velocity perturbations. The rupture time is then directly computed from the rupture velocity distribution, using the algorithm of⁶⁰. This change does not modify the spectral properties of the source functions, but allows setting the maximum rupture velocity to the shear-wave velocity (the prevailing rupture mode is III, which excludes supershear ruptures, e.g., ref. ⁶¹). Figure S3a shows examples of generated rupture models.

**Velocity structure and near-fault wave propagation.** A representative 1D velocity model was derived from the analysis of 6 days of continuous seismic ambient noise recorded synchronously by 19 temporary seismological stations

deployed after the Le Teil earthquake (Fig. 1a, Table S1). Rayleigh and Love wave dispersion curves and signed ellipticity of Rayleigh waves were measured with the 0.1–4 Hz frequency range using a three-component beamforming method (ref. ⁴³) implemented in the open-source Geopsy software (ref. ⁶²) (Fig. S5a-c). Dispersion estimates above 1 Hz are controlled by ambient noise recorded by seismic stations located within 3 km from the epicenter location. Dispersion estimates have been inverted by using the Conditional Neighbourhood Algorithm (ref. ⁴⁴). The best misfit inverted ground velocity profile (Fig. S5d) exhibits a 0.6 km thick layer ($V_s$ ~2 km/s) and a very competent 0.6 km thick layer ($V_s$~3.6 km/s) overlaying a ~1 km thick layer of softer materials. Note that the resolution of the very shallow earth structure (up to 300 m depth) is poor, so that we do not model potential surface rock alteration that may enhance ground acceleration above 5 Hz. Note also that the Triassic bedrock depth is not resolved by the inversion.

The impulse response of the inferred 1D velocity structure (Green’s functions) is next computed for each couple of fault point and virtual station of the dense array using the discrete wave number technique of ⁶³ (AXITA computer package, ref. ⁶⁴). The array is composed of 234 stations with inter-stations distance of 500 m. A distance between the fault points of 60 m ensures stability of the computed ground motion up to 5 Hz.

**Displaced rock analysis.** After the Le Teil Mw 4.9 event, displaced rocks were observed during field surveys on several locations. The fieldwork began 1 month after the earthquake, with several visits within 2 months. A detail description of the observations is provided in Supplementary Information (Supplementary text, Table S2, Figs. S6.1–S6.8). All the observations support the density of upthrown rocks and the location of displaced rocks overlap with the above-gravity peak values of numerical predictions of ground motions (Fig. 3). It supports strong ground motion to be of the order of the gravity during Le Teil earthquake. References on previous worldwide observations of displaced rocks for larger earthquakes are also provided in Supplementary References.

**Observations of funeral slab displacement.** We inspected five cemeteries located in the fault vicinity (Fig. 1a). The structure of French tombs is generally simple, composed of a vault covered by a rectangular granite slab backing onto a headstone (Fig. 5a). We observed that the slabs, which were relying loosely without joints, were displaced during the earthquake. We compile a database of 48 slab displacements from cemetery C1 containing a large majority of the observations and located about 200 m from the northern fault rupture termination (Table S3, Fig. S7). We reported only a few sparse observations with smaller displacements in cemeteries C2 and C3, while no displacement was seen in cemeteries C4 and C5. 75% of the observations in C1 corresponds with slabs backing onto N70°-oriented headstones. Interestingly, within this group of slabs, 80% of the reported displacements are for the slabs backing onto the SE side of the headstones, with an average of 4.5 cm (68% CI: 2–7 cm), against 2.5 cm for the remaining 20% (Figs. S7, S5a).

**Modeling of funeral slab displacement.** We conduct numerical simulations to model the observed funeral slab displacements. We exclude from our analysis slab observations associated with rotation and focus on slabs backing onto N70°-oriented headstones. Since the fundamental frequency of the slab is >100 Hz (considering a granite slab with dimension 200 × 100 × 8 cm), we modeled it as a mass particle with three degrees of freedom. The particle, subjected to gravitational force, is in frictional contact with the EW-NS plane, which is animated by the three-dimensional seismic ground motion. The friction law we used is the standard Coulomb friction model without tangential force regularization. From Newton’s second law, we obtain the differential equation of motion of the undamped particle:

$$ m \left( \sqrt{a_{EW}^2(t) + a_{NS}^2(t)} \right) = F_T(t) \qquad (6) $$

Either $|F_T(t)|<\mu.m(g + a_Z(t))$ if the slab sticks to the EW-NS plane, or $\overrightarrow{F_T(t)} = -\mu.m(g + a_Z(t))$ if the slab moves with a relative velocity $v(t)$ to the plane. $F_T$ is the Coulomb friction force, $\mu$ is the static friction coefficient, $g$ is the constant acceleration gravity, $m$ is the slab mass and $a_{EW}$, $a_{NS}$ and $a_Z$ are ground accelerations in the EW, NS and vertical directions, respectively. During the motion, a perfectly inelastic collision is imposed on the slab, which does not penetrate the headstone (wall azimuth on Fig. 5 and S8, and defined below as azimuth ranges). A forward incremental Lagrange multiplier method compatible with explicit time integration operators (ref. ⁶⁵) is used to solve the particle motion. Convergence is achieved for a time step <10⁻³ s.

We conduct two sets of simulations with two ranges of inhibited azimuth displacement: for set 1 and 2, the slab movement is blocked in the [N-110 – N70] and [N70 – N250] azimuth ranges, respectively. For each set, we consider values of the friction coefficient $\mu$ of 0.2, 0.4 and 0.6 (reported values for granite range between 0.2 and 0.4, ref. ⁶⁶). First, for $\mu = 0.2$, the simulated slab movement is smaller when the azimuth displacement is blocked in the [N-110 – N70] direction (set 1 gives values on average 30% smaller than set 2), as observed on field (Fig. S7). This observation is not reproduced using $\mu = 0.4$ or $\mu = 0.6$ (Fig. S8b). In addition, the values of $\mu = 0.4$ and $\mu = 0.6$ underestimate the displacements observed on field in the SE direction (Fig. S8a). For $\mu = 0.4$ and $\mu = 0.6$, the median observed
displacement (4.5 cm) corresponds to simulated median horizontal PGA values of 9.0 and 12 m/s², respectively, which are not consistent with an Intensity of VIII (ref. ⁵⁵). For μ = 0.2, a displacement of 4.5 cm corresponds to a simulated median horizontal PGA of 4.8 m/s². Accordingly the 0.2 value of friction coefficient μ best approximates the observations. Small variations of μ may explain part of the variability of the field observations.

## Data availability

Seismological data are provided by the French seismological and geodetic network (RESIF Research Infrastructure, http://seismology.resif.fr) including permanent accelerometric data (https://doi.org/10.15778/RESIF.RA), permanent broadband and accelerometric data (https://doi.org/10.15778/RESIF.FR) and seismic noise recorded at temporary post-seismic stations (https://doi.org/10.15778/10.15778/RESIF.3C2019). Other data (databases of displaced objects, velocity model, simulations of ground acceleration and funeral slabs displacements) are available in the Zenodo repository (https://doi.org/10.5281/zenodo.4309183) and in the Supplementary Information (databases of displaced objects, velocity model). The Digital Elevation Models used in Fig. 1a and Fig. 3 are from the Copernicus Emergency Management Service (© European Union, 2012-2020).

## Code availability

All computations and data analyses have been performed using the software Matlab version R2019A (for seismological data analysis, numerical prediction of ground motion, and modeling of slab displacements) and the open-source Geopsy software (ref. ⁶², for noise data analysis and determination of the velocity model). The Matlab codes used to generate individual results are available under request to M.C.

Received: 20 May 2020; Accepted: 18 December 2020;

Published online: 20 January 2021

## References

1. Chiou, B., Darragh, R., Gregor, N. & Silva, W. NGA project strong-motion database. *Earthquake Spectra* **24**, 23–44 (2008).

2. Pacor, F. et al. NESS1: a Worldwide Collection of Strong-Motion Data to Investigate Near-Source Effects. *Seismol. Res. Lett.* **89**, 2299–2313 (2018).

3. Lanzano, G. et al. The pan-European Engineering Strong Motion (ESM) flatfile: compilation criteria and data statistics. *Bull. Earthq. Eng.* **17**, 561–582 (2019).

4. Goulet, C. A., Abrahamson, N. A., Somerville, P. G. & Wooddell, K. E. The SCEC Broadband Platform validation exercise: Methodology for code validation in the context of seismic hazard analyses. *Seismol. Res. Lett.* **86**, 17–26 (2015).

5. Sira C. et al. Rapport macrosismique n°4, Séisme du Teil (Ardèche) 11 novembre 2019 à 11 h 52 locale, Magnitude 5,2 ML (RENASS), Intensité communale max VII-VIII (EMS98), BCSF-RENASS-2020-R2. https://doi.org/10.13140/RG.2.2.27570.84166 (2020).

6. Jomard, H. et al. Transposing an active fault database into a seismic hazard fault model for nuclear facilities – Part 1: Building a database of potentially active faults (BDFA) for metropolitan France. *Nat. Hazards Earth Syst. Sci.* **17**, 1573–1584 (2017).

7. Ritz, J. F. et al. (2020). Surface rupture and shallow fault reactivation during the 2019 Mw 4.9 Le Teil earthquake, France, *Commun. Earth Environ.*, https://doi.org/10.1038/s43247-020-00012-z.

8. Cornou, C. et al. (2020). Rapid response to the Mw 4.9 earthquake of November 11, 2019 in Le Teil, Lower Rhône Valley, *Comptes Rendus Geosci.*, https://doi.org/10.31219/osf.io/3afs5.

9. Klose, D. K. & Seeber, L. Shallow seismicity in stable continental regions. *Seismol. Res. Lett.* **78**, 554–562 (2007).

10. Delouis, B. et al. (2019). Rapport d'évaluation du groupe de travail (GT) CNRS-INSU sur le séisme du Teil du 11 novembre 2019 et ses causes possibles, http://www.cnrs.fr/sites/default/files/press_info/2019-12/Rapport_GT_ Teil_phase1_final_171219_v3.pdf.

11. De Novellis, V. et al. Coincident locations of rupture nucleation during the 2019 Le Teil earthquake. *France and maximum stress change from local cement quarrying. Commun Earth Environ* 1, 20 https://doi.org-insu.bib.cnrs.fr/10.1038/s43247-020-00021-6 (2020).

12. Ampuero, J.-P. et al. The November 11 2019 Le Teil, France M5 earthquake: a triggered event in nuclear country, EGU General Assembly 2020, Online, 4–8 May 2020, EGU2020-18295, https://doi.org/10.5194/egusphere-egu2020-18295, 2020.

13. Pomeroy, P. W., Simpson, D. W. & Sbar, M. L. Earthquakes triggered by surface quarrying-the Wappingers Falls, New York sequence of June, 1974. *Bull. Seismol. Soc. Am.* **66**, 685–700 (1976).

14. Seeber, L., Armbruster, J. G., Kim, W. Y. & Barstow, N. The 1994 Cacoosing Valley earthquakes near Reading, Pennsylvania: a shallow rupture triggered by quarry unloading. *J. Geophys. Res.* **103**, 24505–24521 (1998).

15. Wiejacz, P. & Rudziński, Ł. Seismic Event of January 22, 2010 near Belchatów, Poland. *Acta Geophysica* **58**, 988–994 (2010).

16. Emanov, A. F. et al. Mining-induced seismicity at open pit mines in Kuzbass (Bachatsky earthquake on June 18, 2013). *J. Mining Sci.* **50**, 224–228 (2014).

17. Foulger, G. R., Wilson, M. P., Guayas, J. G., Julian, B. R. & Davis, R. J. Global review of human-induced earthquakes. *Earth-Sci. Rev.* **178**, 438–514 (2018).

18. Grigoli, F. et al. The November 2017 Mw 5.5 Pohang earthquake: a possible case of induced seismicity in South Korea. *Science* **360**, 1003–1006 (2018).

19. Anderson, J. G. Source and Site Characteristics of Earthquakes That Have Caused Exceptional Ground Accelerations and Velocities. *Bull. Seismol. Soc. Am.* **100**, 1–36 (2010).

20. Anderson, J. G. et al. Exceptional Ground Motions Recorded during the 26 April 2008 Mw 5.0 Earthquake in Mogul, Nevada. *Bull. Seismol. Soc. Am.* **99**, 3475–3486 (2009).

21. Bernard, P. & Madariaga, R. A new asymptotic method for the modeling of near-fault accelerograms. *Bull. Seismol. Soc. Am.* **74**, 539–557 (1984).

22. Spudich, P. & Frazer, L. N. Use of ray theory to calculate high-frequency radiation from earthquake sources having spatially variable rupture velocity and stress drop. *Bull. Seismol. Soc. Am.* **74**, 2061–2082 (1984).

23. Schmedes, J. & Archuleta, R. J. Near-Source Ground Motion along Strike-Slip Faults: Insights into Magnitude Saturation of PGV and PGA. *Bull. Seismol. Soc. Am.* **98**, 2278–2290 (2008).

24. Brune, J. N. Tectonic stress and the spectra of shear waves from earthquakes. *J. Geophys. Res.* **75**, 4997–5009 (1970).

25. Kanamori, H. & Anderson, L. Theoretical basis of some empirical relations in seismology. *Bull. Seismol. Soc. Am.* **65**, 1073–1095 (1975).

26. Allmann, B. P., and P. M. Shearer (2009). Global variations of stress drop for moderate to large earthquakes. *J. Geophys. Res.* **114**, https://doi.org/10.1029/2008JB005821.

27. Courboulex, F., Vallée, M., Causse, M. & Chounet, A. Stress-Drop Variability of Shallow Earthquakes Extracted from a Global Database of Source Time Functions. *Seismol. Res. Lett.* **87**, 912–918 (2016).

28. Mueller, C. S. Source pulse enhancement by a deconvolution of an empirical Green's function. *Geophys. Res. Lett.* **12**, 33–36 (1985).

29. Courboulex, F. et al. High-Frequency Directivity Effect for an Mw 4.1 Earthquake, Widely Felt by the Population in Southeastern France. *Bull. Seismol. Soc. Am.* **103**, 3347–3353 (2013).

30. Ross, Z. E., Kanamori, H. & Hauksson, E. Anomalously large complete stress drop during the 2016 Mw 5.2 Borrego Springs earthquake inferred by waveform modeling and near-source aftershock deficit. *Geophys. Res. Lett.* **44**, 5994–6001 (2017).

31. Mordret, A. et al. Seismic stereometry reveals preparatory behavior and source kinematics of intermediate-size earthquakes. *Geophys. Res. Letters* **47**, e2020GL088563 (2020).

32. Heaton, T. H. Evidence for and implications of self-healing pulses of slip in earthquake rupture. *Phys. Earth Planet. Inter.* **64**, 1–20 (1990).

33. Chounet, A., M. Vallée, M. Causse, and F. Courboulex (2017). Global catalog of earthquake rupture velocities shows anticorrelation between stress drop and rupture velocity, *Tectonophysics*, https://doi.org/10.1016/j.tecto.2017.11.005.

34. Tinti, E., Fukuyama, E., Piatanesi, A. & Cocco, M. A Kinematic Source-Time Function Compatible with Earthquake Dynamics. *Bull. Seismol. Soc. Am.* **95**, 1211–1223 (2005).

35. Ohnaka, M. & Yamashita, T. A cohesive zone model for dynamic shear faulting based on experimentally inferred constitutive relation and strong motion source parameters. *J. Geophys. Res.* **94**, 4089–4104 (1989).

36. Gusev, A. A. & Chebrov, D. On Scaling of Earthquake Rise-Time Estimates. *Bull. Seismol. Soc. Am.* **109**, 2741–2745 (2019).

37. Bouchon, M., D. Hatzfeld, J. A. Jackson, and E. Haghshenas (2006). Some insight on why Bam (Iran) was destroyed by an earthquake of relatively moderate size, *Geophys. Res. Lett.* **33**, https://doi.org/10.1029/2006GL025906.

38. Wallace, R. E. Eyewitness account of surface faulting during the earthquake of 28 October 1983, Borah Peak, Idaho. *Bull. Seism. Soc. Am.* **74**, 1091–1094 (1984).

39. Wilkinson, M. W. et al (2017). Near-field fault slip of the 2016 Vettore Mw 6.6 earthquake (Central Italy) measured using low-cost GNSS, *Sci. Rep.* 7, https://doi.org/10.1038/s41598-017-04917-w.

40. McGarr, A. & Fletcher, J. B. Near-Fault Peak Ground Velocity from Earthquake and Laboratory Data. *Bull. Seismol. Soc. Am.* **97**, 1502–1510 (2007).

41. Causse M., F. Cotton and P. M. Mai (2010). Constraining the roughness degree of slip heterogeneity. *J. Geophys. Res.*, https://doi.org/10.1029/2009JB006747.

42. Candela, T., Renard, F., Schmittbuhl, J., Bouchon, M. & Brodsky, E. E. Fault slip distribution and fault roughness: fault slip distribution and fault roughness. *Geophys. J. Int.* **187**, 959–968 (2011).

43. Wathelet, M., Guillier, B., Roux, P., Cornou, C. & Ohnberger, M. Rayleigh wave three-component beamforming: signed ellipticity assessment
from high-resolution frequency-wavenumber processing of ambient vibration arrays. *Geophys. J. Int.* **215**, 507–523 (2018).

44. Wathelet, M. (2008). An improved neighborhood algorithm: parameter conditions and dynamic scaling, *Geophys. Res. Lett.*, **35**, https://doi.org/10.1029/2008GL033256

45. Elmi S. et al (1996) - Notice explicative, Carte géol. France (1/50000), feuille Aubenas (865). Orléans: BRGM, 170 p. Carte géologique par Y. Kerrien (coord.), S. Elmi, R. Busnardo, G. Camus, G. Kieffer, J. Moinereau, A. Weisbrod (1989).

46. Martinez, M. (2013). Calibration astronomique du Valanginien et de l'Hauterivien (Crétacé inférieur): implications paléoclimatiques et paléocéanographiques. Stratigraphie. Université de Bourgogne. Français. tel-00906955v1

47. Aki, K. and P. G. Richards (2002). Quantitative Seismology, 2nd Ed., University Science Books, Sausalito, CA.

48. Aagaard, B. T., Hall, J. F. & Heaton, T. H. Effect of fault dip and slip rate angles on near-source ground motions: why rupture directivity was minimal in the 1999 Chi-Chi, Taiwan, earthquake. *Bull. Seismol. Soc. Am.* **94**, 155–170 (2004).

49. Somerville, P., Smith, N. F., Graves, R. & Abrahamson, N. A. Modification of empirical strong ground motion attenuation relations to include the amplitude and duration effects of rupture directivity. *Seismol. Res. Lett.* **68**, 199–222 (1997).

50. Bouchon, M. et al. Observations of vertical ground accelerations exceeding gravity during the 1997 Umbria-Marche (central Italy) earthquakes. *J. Seismol. 4*, 517–523 (2000).

51. Hough, S. E. et al. Near-Field Ground Motions from the July 2019 Ridgecrest, California, Earthquake Sequence. *Seismol. Res. Lett.* XX, 1–14 (2020).

52. Schlupp A. et al (2020). EMS98 intensity estimation of the shallow Le Teil earthquake, M<sub>L</sub> 5.2, by Macroseismic Response Group GIM. EGU General Assembly, Vienna, Austria, May 3-8, 2020, EGU2020-3767.

53. Caprio, M., Tarigan, B., Worden, C. B., Wiemer, S. & Wald, D. J. Ground Motion to Intensity Conversion Equations (GMICEs): A Global Relationship and Evaluation of Regional Dependency. *Bull. Seismol. Soc. Am.* **105**, 1476–1490 (2015).

54. Bommer, J. J. et al. On the use of logic trees for ground-motion prediction equations in seismic-hazard analysis. *Bull. Seismol. Soc. Am.* **95**, 377–389 (2005).

55. Supino, M., Festa, G. & Zollo, A. A probabilistic method for the estimation of earthquake source parameters from spectral inversion: application to the 2016–2017 Central Italy seismic sequence. *Geophys. J. Int.* **218**, 988–1007 (2019).

56. Abercrombie, R. E., Bannister, S., Ristau, J. & Doser, D. Variability of earthquake stress drop in a subduction setting, the Hikurangi Margin, New Zealand. *Geophys. J. Int.* **208**, 306–320 (2017).

57. Boatwright, J. A dynamic model for far-field accelerations. *Bull. Seism. Soc. Am.* **72**, 1049–1068 (1982).

58. Kaneko, Y., and P. M. Shearer (2015). Variability of seismic source spectra, estimated stress drop, and radiated energy, derived from cohesive-zone models of symmetrical and asymmetrical circular and elliptical ruptures, *J. Geophys. Res.*, **120**, https://doi.org/10.1002/2014JB011642.

59. Dujardin, A. et al (2019). Optimization of a Simulation Code Coupling Extended Source (k⁻²) and Empirical Green's Functions: application to the Case of the Middle Durance Fault. *Pure Appl. Geophys.*, https://doi.org/10.1007/s00024-019-02309-x.

60. Podvin, P. & Lecomte, I. Finite difference computation of traveltimes in very contrasted velocity models: a massively parallel approach and its associated tools. *Geophys. J. Int.* **105**, 271–284 (1991).

61. Freund, L. B. (1990). Dynamic Fracture Mechanics, Cambridge Univ. Press, New York.

62. Wathelet, M. et al. Geopsy: a User-Friendly Open-Source Tool Set for Ambient Vibration Processing. *Seismol. Res. Lett.* **9**, 1878–1889 (2020).

63. Bouchon, M. A simple method to calculate Green's functions for elastic layered media. *Bull. Seismol. Soc. Am.* **71**, 959–971 (1981).

64. Cotton, F. & Coutant, O. Dynamic stress variations due to shear faults in a plane-layered medium. *Geophys. J. Int.* **128**, 676–688 (1997).

65. Baillet, L., V. Linck, S. D'Errico, B. Laulagnet and Y. Berthier. Finite element simulation of dynamic instabilities in frictional sliding contact. In *STLE/ASME 2003 International Joint Tribology Conference* (pp. 25–30). (American Society of Mechanical Engineers Digital Collection, 2003)

66. Biran, O., Hatzor, Y. H. & Ziv, A. Micro-scale roughness effects on the friction coefficient of granite surfaces under varying levels of normal stress. *Meso-Scale Shear Physics in Earthquake and Landslide Mechanics*, (eds. Y. Hatzor, J. Sulem, I. Vardoulakis), (CRC Press, 2009), p. 145–156.

## Acknowledgements

We thank Stéphane Baize, Pierre-Yves Bard, Michel Bouchon, Marc Cushing, Isabelle Douste-Bacqué, Stéphane Guillot, Estelle Hannouz, Sébastien Hok, Christophe Larroque and Jean-François Ritz for fruitful discussions. We thank Romain Jolivet for providing the InSAR slip distribution before publication. This work was partially financed by SINAPS@ research project (ANR-11-RSNR-0022) and by the Ministry of Ecological Transition (research program "Prévention des Risques", N° 2201198104). J.R.G.'s contribution was partially supported by Fonds européen de développement régional (FEDER) (SISM@LP-Swarm research project, POAI-PA0014885).

## Author contributions

M.C. conceived the paper and figures, modeled ground acceleration, analyzed and interpreted farfield seismological recordings and data of modeled slab displacements. C. C. contributed to the acquisition and dissemination of seismic noise data, analyzed noise data, computed and interpreted the velocity model. E.M. designed the field survey in cemeteries, estimated macroseismic Intensity and analyzed databases of accelerometric data. J.R.G. designed the field survey of displaced rocks, compiled the database, and contributed to the discussion section. L.B. modeled the slab displacements. E.E. contributed to the ground acceleration simulations. M.C., C.C., E.M., J.R.G. and L.B. drafted the paper. M.C., C.C., E.M., J.R.G. and E.E. participated to the fieldwork. All authors contributed to the results interpretation and reviewed the paper.

## Competing interests

The authors declare no competing interests.

## Additional information

Supplementary information is available for this paper at https://doi.org/10.1038/s43247-020-00089-0.

Correspondence and requests for materials should be addressed to M.C.

Peer review information Primary handling editor: Joe Aslin

Reprints and permission information is available at http://www.nature.com/reprints

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2021