### DInSAR data analysis

In order to investigate the ground displacements associated with the considered seismic event, we exploited the Differential Synthetic Aperture Radar (SAR) Interferometry (DInSAR) technique, which allows the analysis of surface displacements along the radar line-of-sight (LOS). The SAR data considered were acquired by the Sentinel-1 (S1) constellation of the Copernicus European Program. Benefiting from the short revisit time and the small spatial baseline separation of the S1 constellation, we generated several coseismic differential interferograms (Supplementary Table 2) with a spatial sampling of about 15 m following an averaging operation (multilook) in the azimuth direction to reduce the speckle noise. Among the generated interferograms we selected those less affected by undesired phase artifacts (atmospheric phase delays, decorrelation noise, etc.) for the seismic source modeling discussed in the following section, thus preserving good spatial coverage and interferometric coherence. In particular, the employed S1 data pairs were acquired on November 6–12, 2019 (Supplementary Fig. 2A), and on October 31 and November 12, 2019 (Supplementary Fig. 2B) along the ascending (ASC) and the descending (DESC) orbits, respectively. On both interferograms, several fringes located near the epicentral area are clearly visible; note that each fringe corresponds to a LOS displacement of about 2.8 cm (i.e., half of the employed S1 C-band wavelength *λ* = 5.546 cm). Subsequently, starting from the selected interferograms, we generated their corresponding LOS displacement maps (Fig. 1a) through a phase unwrapping operation^{24}.

By superimposing the fault traces on the selected pair of interferograms (Supplementary Fig. 2C, D and Fig. 1d) it is evident that the deformed area, extending for about 12 km^{2}, is elongated in the NE–SW direction, very similar to the orientation of the LRf. Moreover, a clear asymmetry is clearly visible in the northern sector of the deformed zone, in correspondence of a minor fault located at the base of La Chade Hill’s NW flank (Supplementary Fig. 2E, F), suggesting a complex geometry of the ruptured area. This is also evident in the map of the vertical component of the retrieved deformation pattern (Supplementary Fig. 3A) that we computed from the unwrapped interferograms, together with the East-West deformation map (Supplementary Fig. 3B). Note that the maximum uplift is about 11 cm while the amplitude of the subsidence is appreciably smaller, with a maximum displacement of about 4 cm. For the horizontal displacements, the area affected by westward motion is rather extended, with a maximum value of about 8 cm, while eastward motions are more concentrated on the SW side of the LRf and shows a maximum displacement of about 7 cm.

### Surface deformation data modeling

In order to retrieve the fault parameters, we jointly inverted the DInSAR displacements retrieved from ASC and DESC orbits by applying a consolidated two-step approach that consists of a nonlinear optimization to constrain the fault geometry, assuming a uniform slip, followed by a linear inversion to retrieve the slip distribution on the fault plane^{25}. We modeled the LOS displacements retrieved from the DInSAR interferograms with a finite rectangular dislocation in an elastic and homogeneous half-space^{19}, also applying a compensation for the local topography^{26} and assessing possible offsets and linear ramps affecting the DInSAR measurements. Moreover, data were preliminarily resampled over a regular grid (70 m spacing in all the considered area) to reduce the computation load. Starting from a nonlinear inversion algorithm based on the Levenberg–Marquardt least-squares approach, we searched for the source parameters of one or two rectangular dislocations with uniform slip and, thanks to multiple random restarts implemented within the approach, it was possible to catch the global minimum during the optimization process.

The second step of our modeling is represented by the linear inversion process with the computation of the nonuniform slip distribution, in order to have a more accurate estimate of the slip distribution along the fault plane. In particular, the linear inversion was performed by using as starting model, in terms of dimensions and orientation, the fault obtained from the previous nonlinear inversion discretized into 0.1 × 0.1 km^{2} patches and inverting the following system:

$$left[ {begin{array}{*{20}{c}} {{boldsymbol{d}}_{{boldsymbol{DInSAR}}}} \ 0 end{array}} right] = left[ {begin{array}{*{20}{c}} G \ {k cdot nabla ^2} end{array}} right] cdot {boldsymbol{m}},$$

where *d*_{DInSAR} represents the DInSAR displacements vector, *G* is the Green’s matrix with the point-source functions, ** m** is the vector of slip values for each patch (initially assumed as the value resulting from the nonlinear inversion), and ∇

^{2}is a smoothing Laplacian operator weighted by an empirical coefficient

*k*to guarantee a reliable slip varying across the fault. The choice of the parameter

*k*depends on the compromise between the data fit and the smoothness of the slip distribution. We tested several values and selected

*k*= 0.007, since appreciably higher residuals resulted for

*k*≥ 0.01 and similar residuals but inconsistent slip values (larger than 70 cm), with too rapidly varying slip distribution, resulted for

*k*≤ 0.004. Further constraints were introduced by imposing nonnegative slip (reverse slip only) values obtained via nonlinear inversion.

The surface deformation derived from the interferometric measurements acquired along ASC and DESC orbits reveals a spoon-like geometry, characterized by a NE–SW striking main distinctive displacement pattern (Supplementary Fig. 2). Thus, we first investigated solutions associated with homogeneous dislocation on a single planar source for which all source parameters were set free during the nonlinear inversion. Then, keeping the plane obtained in the geometry inversion fixed, we searched for the best slip distribution. The preferred solution (Supplementary Table 1) consists of a 4.1 × 1.0 km^{2} reverse fault, oriented N50° and dipping 62.3° southeast, with rake 116.5° and dislocation characterized by two separate major patches, located respectively north and south of the quarry, with a maximum slip in excess of 0.3 m (Supplementary Fig. 4). Although the single fault solution accounts for most of the observed DInSAR data, it is associated with large residuals, comparable to the maximum surface deformation (Supplementary Fig. 4). Also, it is worth noting that this solution does not align with any, and actually crosscuts all the faults mapped in the area. Thus, considering these observations, we tested a composite solution, with two fault planes. Similar to what was done for the single fault, we first searched for the best geometry for the two planes, by means of nonlinear inversion with uniform dislocation, and then solved for the slip distribution on the two fault surfaces with fixed geometry by linear inversion. Based on the features of the DInSAR data described above and on the local geology^{11,12}, we tested several positions for the two faults in the uniform slip, nonlinear inversion, always keeping free strike, dip, rake (uniform for each plane), and dislocation. Eventually, we obtained solutions (Supplementary Table 1) noticeably reducing the residuals and, most importantly, compatible with the faults observed at the surface (Fig. 2b and Supplementary Figs. 2 and 3). The two planes coincide respectively with the central sector of LRf and to a structural lineament at the base of the La Chade Hill, both supposed to be SE dipping^{11}, where the quarry is located. The preferred solution is characterized by a slip distribution shallower than 1 km depth. In particular, for the plane *F1* the solution displays a 3 km-long, two-lobed patch with maximum slip of 0.29 m located on the northern half of the fault, while a maximum slip of 0.21 m results for *F2*, located approximately at the center of the fault.

Finally, as for the uncertainties associated with the slip solution we report the standard deviation (Supplementary Fig. 5) as obtained from the model covariance matrix, following an approach^{25} that accounts also for the noise covariance.

### Source directivity analysis and stress drop estimate

We studied source directivity by using a simplified version of the directivity function *C*_{d} (ref. ^{27}) for a bilateral linear rupture model:

$$C_d = frac{1}{2}sqrt {frac{{left( {1 + e} right)^2}}{{left( {1 – alpha rm{cos}vartheta } right)^2}} + frac{{left( {1 – e} right)^2}}{{left( {1 + alpha rm{cos}vartheta } right)^2}}},$$

(1)

where *ϑ* is the angle between the ray leaving the source and the direction of rupture propagation *φ* (ref. ^{28}), and *α* is the Mach-number, that is, the ratio between the rupture velocity *v*_{r} and the S-wave velocity. The percent unilateral rupture *e* parameter is defined as (2*L*′ − *L*)/*L*, where *L* is the total rupture length and *L*′ is the length of the dominant rupture^{29}. A value of *e* = 1 corresponds to a unilateral rupture, whereas *e* = 0 corresponds to a bilateral rupture.

The effect of the directivity on the source spectrum is to increase the corner frequency at those stations located in the same direction as the rupture propagation and to decrease the corner frequency in the opposite direction (e.g., ref. ^{30}). This effect can be modeled assuming that the apparent corner frequency is given by *f*_{ca} = *f*_{c}*C*_{d} with *f*_{c} indicating the actual corner frequency and *f*_{ca} the corner frequencies at various azimuths estimated from the source spectrum.

We retrieved waveforms corrected by the instrumental response at 22 stations from the Réseau Sismologique et Géodesique Français (http://seismology.resif.fr/#WelcomePlace; last accessed April 2020) (Supplementary Fig. 6A). We filtered all the waveforms in the frequency band 0.01–20 Hz. First, we discarded the stations with a low signal-to-noise ratio, which resulted to be poor at all the stations located at an elevation higher than 900 m, reducing the number of stations to 16. Then, we windowed the waveforms by cutting from 2 s before the manual S-wave picking to 6 s after. We applied a 5 per cent cosine taper function and zero padding before computing the Fourier amplitude spectra. The spectra were then smoothed by applying an average moving window with a four-point half width. Finally, we computed the S-wave displacement spectra from the modulus of the three components velocity or acceleration spectra by dividing the spectra by 2*π* or 4*π*^{2}, respectively.

We assumed an *ω*^{−2} theoretical source spectrum model^{31}, which is given by:

$${S}left( f right) = frac{{{{Omega }}_0}}{{1 + left( {frac{f}{{f_c}}} right)^2}},$$

(2)

where Ω_{0} is the long-period spectral amplitude, *f* the frequency, *f*_{c} the corner frequency. As for the anelastic attenuation we assumed the model (Aleft( {f,T} right) = e^{ – pi fT/Qleft( f right)}), where *Q(f)* = *Q*_{0}*f* ^{n} and *T* is the travel time of the *S* phase. Using Eq. (2), together with the attenuation model, we fit the observed spectra through a grid search approach, to infer *Ω*_{0}, *f*_{c}, *Q*_{0} and *n*. As for *Ω*_{0}, we normalized the spectrum and explored the range (0.7, 1.3), while the range of exploration for *f*_{c} was (0.1, 2.0). Finally, we explored *Q*_{0} in the range (0,300) and *n* in the range (0, 2.0). The parameter *Ω*_{0} is then used to compute the seismic moment^{32} through the formulation:

$$M_o = frac{{4pi rho c^3R{mathrm{{Omega} }}_0}}{{FR_{theta varphi }}},$$

where *R* is the hypocentral distance, *ρ* is the density, assumed here 2690 kg m^{−3}; and *c* is S-wave velocity assumed to be 3.5 km s^{−1}. (R_{theta varphi }) is the average S-wave radiation pattern assumed to be 0.7 and *F* is the free-surface coefficient (fixed to 2). As for the uncertainty on each inferred parameter we used an approach^{33} providing confidence intervals based on jack-knife variance analysis. The computed *f*_{c} as function of the station azimuth are finally fit by using the *C*_{d} function and the nonlinear Levenberg–Marquardt least-squares algorithm.

The observed displacement spectra together with the best-fit models are shown in the upper panels of Supplementary Fig. 6B. The average seismic moment value is 3.4 × 10^{16} (2.4 × 10^{16}, 4.9 × 10^{16}) N m corresponding to *M*_{W} = 4.9. As for the anelastic attenuation, we found that a *Q* frequency independent model (i.e., *n* = 0) provides the best fit, with *Q*_{0} = 190 (193, 198). The estimated corner frequencies are shown in the lower panel of Supplementary Fig. 6B as a function of the station azimuth. The result clearly indicates at least one range of azimuths where the *f*_{c} value increases with respect to the other directions. The fit with the *C*_{d} function, which is shown in the same panel, provides an *e*-value of 0.3 ± 0.1 indicating a bilateral rupture and a Mach-number of 0.5 ± 0.1 suggesting a rather low rupture propagation. The dominant rupture direction is at 241 ± 8° while the secondary direction is at about 60°. Using the average value of the estimated corner frequencies corresponding to 0.6 (0.5, 0.7) Hz and the seismic moment, we computed the static stress Δ*σ* = 0.44 *M*_{0}/*r*^{3}, with the source radius given by *r* = 0.37*v*_{S}/*f*_{c} (ref. ^{28}), being *v*_{S} the S-wave velocity. We obtained *r* = 2158 (2044, 2846) m and Δ*σ* = 1.3 (0.9, 2.0) MPa.

### Calculation of the rock volume extracted from the quarry

We used multi-temporal DSM to estimate the removed volume of rock in the Le Teil quarry. Topographic map resolution, quality, and accuracy varies a lot with time, so a more robust method of change detection using high-resolution topography from stereo aerial imagery is preferred. However, the Le Teil quarry was established in 1833 and this technique is only available for the modern era. Thus, distinct methods must be used for the different time periods.

For the recent period, we used scanned stereo aerial images of the study area, available from the Institute National de l’Information Géographique et Forestiére (IGN). We used archive stereo aerial image pairs for 1946, 1979, 2007, and 2011 (Supplementary Table 3). We extracted DSM from stereo imagery using MicMac software^{34}, following the procedure illustrated in Supplementary Fig. 7. For each group of stereo images, tie-points between images and camera frame/lens parameters were calculated. Using ground control points, we performed bundle adjustment and refinement of the camera calibration. Ground control points were picked from the latest IGN orthophoto map (coordinates in UTM zone 31 north, WGS84) and elevation values were extracted from ALOS Global Digital Surface Model. Then, DSM and orthophotos were produced. A point cloud was extracted from the Malt DSM and colored with RGB values from the orthophoto map (Supplementary Fig. 7). A detailed photogrammetric point cloud of 2007 covering a much larger area is available from the OpenTopography facility (https://doi.org/10.5069/G9BC3WQ4; last accessed April 2020).

Incidentally, we remark that there is a strip of images dated 1932 available from IGN, but processing produced many artifacts due to scanned film’s poor preservation and noise. Coregistration with 1946 data was also poor, with large residuals; thus, we decided not to use the 1932 dataset. We did not use the most recent set of digital aerial images provided by IGN (dated 2013), as the poor radiometric quality of the files leads to a low quality matching over the quarry area, thus rendering the dataset unusable for terrain change detection.

In order to estimate the initial reference topography of the site as best as possible, we searched for the earliest available topographic map on a suitable scale. We used the Carte de l’État-Major-feuile Privas S.O. (Supplementary Fig. 8), a map compiled between 1846 and 1857 at a scale 1:40,000. We georeferenced a high-resolution scan obtained from IGN (https://www.geoportail.gouv.fr/donnees/carte-de-letat-major-1820-1866; last accessed April 2020). A concise description of the main map-series features and how to spatially process each map sheet is also available^{35}. As the map does not provide detailed and accurate contour lines, we extracted a dense set of elevation contours from the 1946 DSM. These contours were clipped at the maximum extent of quarry boundaries traced from the 1946 orthophoto (Supplementary Figs. 8 and 9). Using the 1846–1857 map as a reference for general relief representation and taking into consideration the geomorphologic features (Rhône river escarpment, drainage, and ridgelines) of the area around the quarry, contours were traced manually in order to fill the blank area and reconstruct the 1833 relief. By assuming for the reconstructed topographic height an uncertainty of 10 m over the whole quarry area as of 1946 (~300,000 m^{2}), the uncertainty in the volume removed in the period 1833–1946 (larger than 11,000,000 m^{3}) deriving from the procedure of reconstruction of the initial topography would be <30%.

Starting from the evolution of the topography, we estimated the volume of the mass progressively removed from the quarry. Instead of comparing the first and last DSM available, in order to increase coherency and avoid very large outliers we split the analysis into intervals: 1833–1946, 1946–1979, 1979–2007, and 2007–2011. Original DSM raster files were clipped and resampled into a common pixel size (3 m for 1833–1946 and 2 m for the rest) and then converted into point clouds. Point clouds were co-registered using the ICP matching tool from CloudCompare software (http://www.danielgm.net/cc; last accessed April 2020) and then each pair was processed with the Multiscale Model to Model Cloud Comparison (M3C2). This technique computes the local distance between two-point clouds along the normal surface direction which tracks 3D variations in surface orientation^{36}. M3C2 distance raster files were cropped using the quarry boundary traced manually from orthophoto maps for each period, and removed volume (Supplementary Table 4. The uncertainty in the determination of the removed volume for each period is also reported) was calculated using raster statistics.

### Boussinesq’s solution of 3-D elastostatic loading

We adopted the analytic solution in Cartesian coordinates (*x*_{1}, *x*_{2}, *x*_{3}) of the three-dimensional elastostatics for the response of a semi-infinite solid (*x*_{3} ≥ 0) to an arbitrary normal load *P*_{0} on its boundary^{37}. We considered the displacements and the stress components given the shear modulus *μ* and the Poisson ratio ν. If *P*_{0} is applied at (*x*_{01}, *x*_{02}*, x*_{3}) then the displacement components in a generic point (*x*_{1}, *x*_{2}, *x*_{3}) are given by:

$$u_1 = frac{{P_0}}{{4pi mu }}frac{{left( {x_1 – x_{0_1}} right)}}{R}left[ {frac{{x_3}}{{R^2}} – left( {1 – 2nu } right)frac{1}{{R + x_3}}} right],$$

$$u_2 = frac{{P_0}}{{4pi mu }}frac{{left( {x_1 – x_{0_2}} right)}}{R}left[ {frac{{x_3}}{{R^2}} – left( {1 – 2nu } right)frac{1}{{R + x_3}}} right],$$

$$u_3 = frac{{P_0}}{{4pi mu }}frac{1}{R}left[ {frac{{x_3^2}}{{R^2}} + 2left( {1 – nu } right)} right],$$

where (R = sqrt {left( {x_1 – x_{01}} right)^2, +, left( {x_2 – x_{02}} right)^2, +, x_3^2}). The components of the displacement can be used to compute the strain tensor *ε* and the stress tensor σ in the same point (*x*_{1}, *x*_{2}, *x*_{3}). Under the hypothesis of infinitesimal deformation and for isotropic medium, *ε*_{ij} = *1/2(dU*_{i}*/dx*_{j} + *dU*_{j}*/dx*_{i}*)* and *σ*_{ij} = *λ*δ_{ij}*ε*_{kk} + 2*με*_{ij}, where *λ* is the second Lame’s parameter and δ_{ij} the Kronecker delta. In particular, for the application presented in this study we used the six components of the stress tensor, which have the following expressions:

$$sigma _{11} = – frac{{P_o}}{{2pi R^2}}left{ {frac{{3left( {x_1 – x_{01}} right)^2x_3}}{{R^3}} – left( {1 – 2upsilon } right)left[ {frac{{x_3}}{R} – frac{R}{{R + x_3}} + frac{{left( {x_1 – x_{01}} right)^2left( {2R + x_3} right)}}{{Rleft( {R + x_3} right)^2}}} right]} right},$$

$$sigma _{22} = – frac{{P_o}}{{2pi R^2}}left{ {frac{{3left( {x_2 – x_{02}} right)^2x_3}}{{R^3}} – left( {1 – 2upsilon } right)left[ {frac{{x_3}}{R} – frac{R}{{R + x_3}} + frac{{left( {x_2 – x_{02}} right)^2left( {2R + x_3} right)}}{{Rleft( {R + x_3} right)^2}}} right]} right},$$

$$sigma _{12} = – frac{{P_o}}{{2pi R^2}}left( {x_1 – x_{o1}} right)left( {x_2 – x_{o2}} right)left[ {frac{{3x_3}}{{R^3}} – left( {1 – 2upsilon } right)frac{{left( {2R + x_3} right)}}{{Rleft( {R + x_3} right)^2}}} right],$$

$$sigma _{13} = – frac{{3P_o}}{{2pi R^5}}left( {x_1 – x_{01}} right)x_3^2,$$

$$sigma _{23} = – frac{{3P_o}}{{2pi R^5}}left( {x_2 – x_{02}} right)x_3^2,$$

$$sigma _{33} = – frac{{3P_o}}{{2pi R^5}}x_3^2.$$

In order to quantify the effect of the removal of the rock mass from the quarry on the prescribed fault, we assumed that extraction operation is equivalent to the application of a set of vertical forces *P*_{0} = *h*·*dA*·*ρ*·*g* (being *h* the height of the eroded rock column computed from the DEM, *ρ* the rock density, *g* the acceleration of gravity and *dA* the elementary area of the DEM), at the surface of the investigated volume. Given the linearity of the model, the net effect on the fault in terms of *σ*_{ij} is obtained by summing the contributions of the single forces. By decomposing *σ*_{ij} in its normal component Δ*σ*_{n} (assumed positive if the fault is unclamped) and shear component Δ*τ* (assumed positive in the direction of the slip) along the fault and assuming a given value of the fault effective friction coefficient *μ*′ we computed the Coulomb stress change ΔCFF = Δ*τ* + *μ*′Δ*σ*_{n}. Specifically, we investigated three values of *μ’* = 0.4, 0.5, 0.6 (Supplementary Fig. 10), generally assumed as plausible effective friction coefficient for natural faults^{38,39} and minor differences resulted in the maximum values, while no change in the general ΔCFF pattern. Conservatively, in order to estimate the effect of the quarry activity on the faults associated with the Le Teil earthquake, we used the results for *μ* = 0.4 (Fig. 3c), associated with the lower maximum value.