Development of an empirical ground-motion model for post-Development of an empirical ground-motion model for post-mining induced seismicity near Gardanne, France mining induced seismicity near Gardanne, France

Abstract Since the closure of mining activities in 2003, the coal basin of Gardanne in South-East France has experienced thousands of small-magnitude earthquake events, mostly triggered by the flooding of mine workings. Some of these events have been powerful enough to be strongly felt by the population, generating nuisance and concern about potential damage to buildings. The aim of this study is to improve the characterisation of the level of ground motion at the surface, by developing a ground-motion model for post-mining induced seismicity, based on several years of recorded data. A Bayesian-based method is applied to the data in order to account for uncertainties in the estimation of moment magnitude. Station-to-station ground-motion site terms are also quantified for the nine recording stations in the area, thus providing additional information on the local site conditions. The developed model is compared to existing prediction equations for seismicity induced by other types of anthropic activities, confirming the need for a specific model in the case of post-mining induced seismicity. Finally, the Gardanne ground-motion model is also integrated with a shake-map procedure, showing how this predictive model may be merged with recorded data in order to generate rapid estimates of shaking levels in the area.


Introduction
T he management of post-mining areas in Europe is an important issue for safety and economic reasons.Within the PostMinQuake European RFCS research project (https:// postminquake.eu/),potential environmental threats due to mine closures are investigated via the impact of mine flooding on underground mining excavations and stability of rock masses, and its potential to trigger superficial earthquakes [1,2].Ground vibrations caused by seismic events as well as continuous and discontinuous surface deformations (subsidence and sinkholes) can influence the evolution of underground excavations, ground stability and surface infrastructure.
In South-East France, the Gardanne coal basin ceased its mining activities in 2003, and repeated flooding of the mining galleries due to fluctuating levels of the aquifer have been triggering earthquakes, with alternating periods of low and high activity.It is therefore necessary to better characterise the hazard levels to which the population may be exposed, by improving the prediction of the distribution of ground-motion parameters at the surface.
Such earthquakes are usually extremely shallow (less than 1 km) and of very low magnitudes [3], so that the affected area seldom exceeds several square kilometres at the surface.These specificities prevent the application of current predictive models in seismic hazard and risk analyses, which are mostly focused on natural seismicity.In the case of classical ground-motion models (GMMs), which estimate the distribution of ground-motion parameters given some earthquake characteristics, the range of source-to-site distances, magnitudes and rupture mechanisms are not compatible with the seismic events induced by post-mining sequences.For instance, the GMM by Atkinson [4], which is referred to as a GMM adapted to induced seismicity [5], is based on smaller magnitude events (M w between 3 and 6) selected from the NGA-West2 database of tectonic earthquakes [6].However, it should be noted that the lower magnitude limit of this GMM is still way above the range of magnitudes that may be considered in the case of post-mining induced events.More recently, some GMMs have been specifically developed in the context of induced seismicity, however they are mostly focused on earthquakes triggered by geothermal activities [7], hydraulic-fracture operations [8] or natural gas extraction [9].Their applicability to the post-mining induced seismicity of the Gardanne coal basin remains dubious and should be investigated.
Therefore, the aim of this study is to derive an empirical GMM adapted to the Gardanne area, by taking advantage of several years of recordings of post-mining induced events.In this coal basin, Ineris (French National Institute for Industrial Environment and Risks) operates, on behalf of BRGM (French Geological Survey), a permanent seismic network for characterising the recorded events and publishing information bulletins [10].The aim of this survey is to prevent mining mechanical instabilities associated with room-andpillar sectors below inhabited areas, in mining hazard zones qualified by GEODERIS [11,12], which is a Public Interest Group expert in post-mining issues.This permanent seismic network was designed to detect and monitor the first signs of instability at the level of the mining structures and to anticipate potential disorders on the surface.
Since 2010, the basin has been periodically affected by seismic activity unexpectedly located outside of the identified areas at risk of mechanical instabilities, monitored by the permanent seismic network.This seismicity was regularly felt by the population.To better understand the origin of this seismicity, and following a public request, BRGM also operates a network of 9 temporary stations in the immediate vicinity of most of the events, as described in Section 2. Ineris subsequently supplemented this network with 5 additional temporary seismic stations.Then, the first challenge is to consolidate a harmonized catalogue of events and records, in terms of moment magnitude and epicentre location and depth, between the Ineris and BRGM respective datasets.To this end, it is proposed to back-calculate some events with the SourceSpec code [13] and to check potential uncertainties in the moment magnitude estimates (see Section 3).For the derivation of the GMM in Section 4, the Bayesian-based regression method introduced by Kuehn & Abrahamson [14] is applied here, because it allows for the inclusion of uncertainty on predictor variables (e.g., moment magnitude) and for the computation of so-called station-to-station ground-motion site terms [15], which may be seen as proxies for soil amplification effects in the absence of high-resolution site data.Finally, in Section 5, it is shown how the derived GMM can be integrated into a shake-map procedure (i.e., estimation of the ground-motion field by accounting for recorded ground-motion parameters at nearby stations), in order to improve the prediction of the effects of the induced earthquakes and to manage post-mining risks better.

Description of the Gardanne area
The coal basin of Gardanne is located about 20 km NNE of Marseille (Bouches-du-Rhône, France) and it covers ~60 km 2 (see Fig. 1).The industrial exploitation has lasted from the first half of the 19th century to 2003.With more than 500 km of tunnels and levels exploited up to 1350 m deep, the Gardanne mine is the largest mining area in South-East France.During the exploitation, the miners had to pump the underground water.After 2003, the aquifer recovered rapidly until 2010, when new pumps were installed to stabilize the level of underground water between À30 m ASL (above sea level) and þ10 m ASL (at ~300 m depth).In order to prevent a potential hazard of rupture of pillars (which may generate small earthquakes), accompanied by surface effects (sinkhole, subsidence), Ineris started operating a microseismic monitoring permanent network of 5 stations (see Fig. 2 left) on behalf of BRGM.
After a seismic sequence strongly felt at the end of 2012, outside of the areas monitored by the permanent network, 4 new seismic monitoring stations were deployed by BRGM in 2013.As of 2018, the network has been densified (now 9 BRGM stations over around 4 km 2 , as described in Table 1).The microseismicity analysis on the December 2014 sequence suggests an activation of an adjacent, deeper fault system [3,16e18].From September 2016 to April 2017, the activity increased sharply, prompting many reactions from the population of the Fuveau and Gr easque areas.The activity has then decreased, but it persists with a dozen seismic events per month (M À0.5 to 1.5; depths < 1 km) [3] and around 5 events felt each year.Five complementary stations were set by Ineris at the end of 2018 and at the beginning of 2019 in the same area, but they have not been used in this study.Only BRGM stations are used in the ground-motion dataset, in order to preserve homogeneous data in terms of instrumental response and station set-up (i.e., Ineris velocity sensors are located at 2 m depth in a borehole, while BRGM acceleration sensors are at surface level).
Since the deployment of the BRGM seismic network, the highest Peak Ground Acceleration (PGA) has been recorded by station ROSS, at 74 mg, during the October 12th 2016 event (local magnitude in BRGM catalog M l (BRGM) ¼ 1.5).
Over the period 2018e2022, for which the recordings from all the 9 BRGM stations are available  (see Table 1), the maximum PGA has been measured by station BULL on the April 19th 2019 event (M w,Ineris ¼ 1.7; M l,BRGM ¼ 1.3), with a value of 55 mg.For the same event and same station, Peak Ground Velocity (PGV) reaches 0.49 cm/s, which is also the maximum over the 2018e2022 period.Most PGV/ PGA ratios are very low (i.e., around 0.005e0.015s), which implies a majority of high-frequency signals.This observation is in line with the small-magnitude events recorded at very short distances [19].As an example, some recorded spectral curves for Fourier amplitude, fitted with a Brune model [20], are shown in the Appendix (Fig. A1,c).BRGM stations have a sampling rate of 1000 Hz, with a corner frequency of 200 Hz for the sensor and 500 Hz for the digitizer, which is sufficient to record high-frequency signals and to measure the true PGA.

Characterisation of the recorded events: magnitude and epicentre location
Ineris establishes an operational catalogue of the seismicity, and then sends an information bulletin in case of potentially felt earthquakes (duration magnitude greater than 1.8).This constitutes the operational event catalogue, where the moment magnitude and epicentre location are estimated using the Ineris seismic permanent network (see Fig. 2

left).
In order to check the consistency between the procedures used at BRGM and Ineris, 20 earthquakes are selected from the Ineris operational catalogue between January 2019 and July 2019, which are also well recorded by the BRGM seismic network (Fig. 2 right).All the earthquakes are located within the network, except for one event (2019/01/10 16:53).The range of local magnitude is estimated between M l ðBRGMÞ ¼ À0.26 and 1.35 in BRGM catalogue [3], and M D ðInerisÞ ¼ 0.9 and 2.3 in Ineris catalogue.
The 20 events are analysed with the SourceSpec code [13,21], which is based on the Fourier amplitude spectrum of the S-wave displacement in the far field, including geometric and anelastic attenuation of body waves [20,22].The Fourier amplitude spectrum Sðf Þ averaged over the components is fitted by: Where GðrÞ is the geometrical spreading coefficient [23] (functional form is taken as G(r) ¼ r n with n ¼ 1 as default value [13]), R QF is the radiation pattern coefficient (averaged here without considering focal mechanism ¼ 0.62 [13]), r h and r r are the medium densities at the hypocenter and at the station, M 0 is the seismic moment, f c is the corner frequency, t * is the attenuation parameter.The moment magnitude is then given as M w ¼ 2 3 ðlog 10 ðM 0 Þ À 9:1Þ [24].The code estimates simultaneously the three parameters M 0 ðM w Þ; f c and t * .Furthermore, according to the formulation of Madariaga [25], the apparent stress drop Dt and a source radius r are estimated.The quality factors of attenuation Q P and Q S can be estimated as Q P ¼ tt P ðrÞ t * and Q S ¼ tt S ðrÞ t * for P-and Swaves, respectively.The variables tt p and tt s are the travel times of P-and S-waves, respectively.
The data are stored in SAC format with event information (origin time and hypocenter location).A homogeneous elastic medium is assumed, with the following characteristics: density r ¼ 2439 kg/m 3 , Pand S-wave velocities as c h ¼ 3.85 km/s and c r ¼ 2.03 km/s (e.g., [3]).These are used for determining the constant of the above equation and also for estimating the synthetic S-wave arrival at each station.The window is then taken from 0.1 s earlier than the estimated S-wave arrival time for a  In order to explore the statistical features of the estimated parameters among the 20 earthquakes, the different parameters determined in the Sour-ceSpec code are compared (Fig. 3).In the code, the moment magnitude is well estimated when compared to the local magnitude that was independently obtained.The two other parameters f c and t * are compared: although it might have been expected that the corner frequency f c would increase for smaller magnitudes, no clear dependency on magnitude is observed.The attenuation parameter t * seems to decrease for a smaller magnitude in panel (c), but it is also related to the fact that the available station distance is shorter for smaller magnitude as seen in panel (d ).At such short distances, in the 1e3 km range, t* becomes equivalent to k 0 (i.e., the intercept in distance-dependent k models), with values ranging from 0.01 to 0.025 s: these values are somewhat in line with the k model that was estimated by Douglas et al. [26] for various regions in France (e.g., k 0 around 0.02 s for rock sites in South-East France).The quality factor of attenuation is then roughly estimated as: The stress drop (Dt) and a source radius are also estimated, and they might have been strongly influenced by the previous estimation of f c and t * .Expected observations might have been no magnitude dependency in Dt and a positive scaling between source radius and magnitude.Results from the limited dataset apparently show a positive relation with magnitude.The magnitude dependency of stress drop has been debated for a long time (e.g., Abercrombie & Leary [27]; Ide & Beroza [28]).It seems that the positive dependency is often found in a limited dataset; however, this tendency may disappear when compiling different datasets and taking account of missing areas [28].In general, small earthquakes with a high stress drop lead to a high corner frequency, but this is constrained by the observational limit.Moreover, actual observations imply that one can reasonably estimate the amplitude level of spectral, while it is difficult to estimate the corner frequency f c (see Fig. A1).
Then, the obtained moment magnitude M w is compared with the reported magnitudes in the Ineris catalogue (Fig. 4).The results are fitted with a linear equation only for the moment magnitude of the Ineris catalogue equal to or larger than 0.3 (M w ðInerisÞ !0:3), as data are selected with this criterion in the following analysis of GMM.The M w shows a good correlation with local magnitudes independently obtained (panels a and b), regardless of the trend.The obtained M w is found between the two magnitudes reported in different catalogues (panel b and c) with the intercept B ¼ 0.480 and B ¼ À0.689, respectively.A linear relation may be fitted between the two independently obtained moment magnitudes (panel d ).However the direct difference may be more important than the fitting coefficients, as the sampling number is limited.Defining D ¼ M w ðInerisÞ À M w ðSource SpecÞ, the average D is À0.04 with a standard deviation of 0.17.The maximum jDj is 0.42 for the event of 2019/01/16 23:09.Figure 5 illustrates the obtained moment magnitude and the difference D on a map of the 20 events.No clear tendency of magnitude difference in geography can be observed.The event of 2019/01/ 16 is located around 200 m south of one station of Ineris.In the SourceSpec analysis, this earthquake was recorded by all the nearby stations.It is possible that the near-field effect of earthquake source (radiation) affected the estimation.
Finally, the result of this analysis can justify the choice of the full GMM dataset (see Section 4) with the moment magnitude that has been estimated by Ineris: in order to focus only on significant events, only events with moment magnitude M w (Ineris) !0.3 are selected.Moreover, the location of the events has been back-calculated by BRGM [3], showing a potential deviation of up to several hundred of meters between the locations estimated by BRGM and Ineris.This discrepancy is due to the BRGM network being denser and closer to most earthquake locations in this area, compared to the Ineris one.A sensitivity study has been carried out in order to estimate the impact of a location deviation on the estimation of the moment magnitude, showing that a deviation of around 500 m may result in a magnitude error that is normally distributed within a variation of 0.1 (Fig. 5c).In order to limit this effect, seismic events with a location deviation above 1000 m are removed from the GMM dataset.

Derivation of the ground-motion model
An empirical GMM results from the regression on recorded data points, i.e. the estimation of a groundmotion parameter as a function of predictor variables, such as the event magnitude or the source-to-site distance.Due to the scarcity of existing GMMs for induced seismicity [29], and especially for seismicity following coal mine closures; this section focuses on the development of a specific GMM for the Gardanne site by exploiting the data collected since the end of mining activities.

Catalogue of ground-motion data
The preliminary work on the recorded events (see Section 3) has led to the selection of 94 seismic events recorded by 9 stations operated by BRGM (see Table 1), for the period 2018e2022, as shown in Figure 6.
As a result, 539 acceleration time-histories recorded by the 9 stations have been subjected to the following processing steps: Baseline correction and integration in order to extract the velocity time-history.Estimation of a signal-to-noise ratio for each time-history (i.e., computation of a filtered envelope of the time-history, and evaluation of the ratio as the maximum envelope value divided by the first value of the envelope) and removal of the low-quality recordings.Extraction of the ground-motion parameters PGA and PGV, designated as the geometric mean of the peak values of the two horizontal components of the acceleration and velocity timehistories, respectively.Computation of the spectral acceleration (SA) at periods T ¼ [0.02 s; 0.05 s; 0.1 s; 0.2 s; 0.3 s; 0.5 s], taken as the geometric mean of SA of the two horizontal components of the time-histories.Due to the very high-frequency content and short duration of these specific signals, it has been found that the acceleration spectrum at larger periods was not exploitable, preventing the derivation of a GMM for SA(1.0 s), for instance.
The distribution of the 539 data points, in terms of moment magnitude and hypocentral distance (defined as the distance to the earthquakes located by BRGM), is plotted in Figure 7. Due to the very specific seismic context, a limited range of magnitude (M w between 0.3 and 1.7) and hypocentral distance (R hyp up to 7.5 km only) is available.

Construction of empirical GMMs
For a ground-motion parameter of interest Y (e.g., PGA, PGV or SA), the GMM usually takes the following form [14]: where Q are the coefficients of the model, x is a vector comprising the predictor variables, and the error terms are: dB, the between-event residual, following a normal distribution of standard deviation t; dS, the site-to-site residual, following a normal distribution of standard deviation f S2S ; dWS, the within-event and within-site residual, following a normal distribution of standard deviation f SS .
The functional form f(Q,x) of the GMM has been selected from a series of statistical trials on several functional forms, inspired by existing references such as Joyner & Boore [30], Douglas et al. [7] or Atkinson [4].A quadratic magnitude scaling term has been added to the model in order to better constrain predicted values at larger magnitudes and to avoid unrealistic results if the GMM is used for magnitudes that are slightly above the range provided by the dataset (i.e., M w up to 1.7).Regarding distance-dependent terms, no anelastic decay term  is included due to the absence of observations outside the near-source zone.Although the chosen distance metric is hypocentral distance, it has been decided to introduce a term log 10 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R 2 hyp þ h 2 q in the GMM, with h ¼ 0.1 km as a fixed constraint.This proves to be a satisfying compromise between the compliance to empirical data and the need to saturate predicted values in the rare cases where R hyp is less than 0.1 km (the minimum depth of recorded events is around 150 m).Therefore, the final functional form is the following, including magnitude dependence of distance scaling: where Q ¼ [c 1 ; c 2 ; c 3 ; c 4 ; c 5 ] and x ¼ [M w ; R hyp ] are respectively the vectors of GMM coefficients and predictor variables.The parameter h is set at 0.1 km.Due to the uncertainty of the magnitude estimation, the Bayesian-based method by Kuehn & Abrahamson [14] for the derivation of the GMM is used here.It allows to take into account measurement uncertainty in predictor variables, such as magnitude, in the development of a GMM.The model treats the value of the predictor variable as a parameter to be estimated, constrained by its observed value and the estimated variance of the measurement error.The observed (or estimated) moment magnitude is assumed to follow a normal distribution, where the mean M w * is the unknown true value and the standard deviation is s Mw : The standard deviation s Mw is set to 0.1, based on the aforementioned sensitivity study on the influence of the location error on the magnitude estimation.
In the Bayesian regression applied by Kuehn & Abrahamson [14], the objective is to estimate the posterior distribution of the predictor variables Q, as well as of the values of the standard deviations t, f SS , and f S2S , and the values of the event terms dB and stations terms dS for each event and station.
Given the set of parameters L to estimate and D the recorded dataset, the posterior distribution of L is expressed as follows, according to Bayes' rule: where P(DjL) is the likelihood of the data and P(L) is the prior distribution.
The prior distributions of the predictor variables are defined as informative prior distributions, as recommended by Kuehn & Abrahamson [14].Their values are taken from the results of a preliminary GMM regression (i.e., non-linear regression based on a gradient search of the GMM coefficients that match the recorded dataset).
Finally, the Bayesian updating is implemented and performed in the Stan code [31] by following the modelling template shared by Kuehn & Abrahamson [14].The posterior distributions of the parameters are obtained by using Markov Chain Monte Carlo sampling: 4 chains with 10 000 samples for each chain are run, with an initial burn-in period of 2000 samples.The sampled values are then used to build an empirical posterior distribution of the GMM parameters.The Bayesian regression each time separately for each of the studied groundmotion parameters.The regression results are detailed in Table 2.
The derived GMMs for Gardanne are plotted in Figure 8 (PGA) and Figure 9 (PGV) for some earthquake scenarios (M w 0.5 and 1.0).They are compared to the data points that have been used in the regression (selected within a À/þ 20% interval around the magnitude of interest).
In case of the PGA model, the Gardanne GMM is compared to three existing ones, which have been selected from the GMM compendium of Douglas [29]: McGarr & Fletcher [32] GMM: model developed from around one hundred of records generated  M w between 1.0 and 4.0, depth up to 10 km).Cremen et al. [8] GMM: use of records from earthquakes induced by hydraulic-fracture operations, due to shale gas development in the UK (hypocentral distance between 2 and 7 km, M w between 0.0 and 3.0, depth between 1000 and 2500 m).While these GMMs have not been specifically developed for post-mining induced seismicity, they cover similar types of events in terms of magnitude and distance.Regarding focal depth, the GMM by McGarr & Fletcher [32] is based on similar ranges as the Gardanne events, while the ones from Douglas et al. [7] and Cremen et al. [8] are based on somewhat deeper events.It appears that the GMM by McGarr & Fletcher [32] underestimates the PGA predictions by around half an order of magnitude, which justifies the development of an ad hoc GMM specifically for the Gardanne area.The GMM by Cremen et al. [8] has a very similar functional form, but it slightly overestimates the PGA predictions and it has a larger decay rate than the Gardanne GMM: this overestimation may be due to the fact that the Cremen et al. [8] model has been derived for alluvium sites (i.e., V s,30 around 280 m/s), while the Gardanne GMM is mostly designed for stiff soil or rock conditions.Regarding the Douglas et al. [7] GMM, it appears that the applied near-source distance saturation term is not adapted to the Gardanne dataset, leading to a substantial underestimation of PGA predictions at distances less than 2 km.This much flatter slope observed for the Douglas et al. [7] GMM is likely due to the larger depth of events used by their model (i.e., up to 10 km depth) when compared to the shallow depths of the recorded events in the Gardanne basin.However, it should be noted that, within the 2e4 km range of hypocentral distances, the Gardanne GMM and the models by Douglas et al. [7] and Cremen et al. [8] still provide comparable estimates.In terms of variability, the total standard deviation of the proposed GMM (s tot ¼ 0.366) matches the values from similar models, such as Atkinson [4] (s tot ¼ 0.370), McGarr & Fletcher [32] (s tot ¼ 0.483) or Cremen [8] (s tot ¼ 0.309).Fig. 8. GMM for PGA (red lines) with two magnitudes (M w 0.5 and M w 1.0), compared with the data points used in the regression, and with a GMM from the literature.The dotted lines represent the þ/À 1 standard deviation boundaries.
Regarding the inter-event variability, the standard deviation of the proposed GMM (t ¼ 0.291) is also comparable to the ones found in the literature, e.g.t ¼ 0.240 in Atkinson [4] or t ¼ 0.190 in Cremen et al. [8].The variability in the Douglas et al. [7] GMM is much larger, which may be due to the wide range of geographical areas and seismicity contexts that were covered in their dataset.
In the case of PGV predictions (see Fig. 9), the Gardanne GMM is only compared to Douglas et al. [7] and Cremen et al. [8] GMMs, since the model from McGarr & Fletcher does not provide PGV estimates.Globally, the same observations as in the PGA case may be reached, with the difference that the decay rate from the Cremen et al. [8] GMM is even higher relatively to the one from the Gardanne GMM.
The spectral shape obtained from the GMM is also investigated in Figure 10.Two specific scenarios have been chosen (M w 1.0 and R hyp ¼ 1 km; M w 0.5 and R hyp ¼ 1.5 km) due to the relatively large number of data points available at these magnitude/ distance combinations.This selection allows to compare the recorded spectral shapes with the ones predicted by the GMM (Fig. 10, left and middle), showing a very good agreement both for the mean and for the standard deviation.In Figure 10 right, the period slope of the normalized response spectrum is compared for various magnitudes (M w 0.0, M w 1.0, and M w 1.5).For longer periods well below the corner frequency, the slope should be 2 with an omega-squared point-source model.Here, the slope is found to be very close to 2 for the considered magnitude range, thus demonstrating the consistency of the developed GMM across various periods.It is also observed that, for M w 0.0 and M w 1.0, the short-period spectral shape appears to be influenced by the attenuation parameter t*: this may be due to the low values of t* at such magnitudes (see Fig. 3,c) leading to an attenuation corner frequency f k ¼ 1/(pt*) around 30 Hz, which is slightly superior to the estimated corner frequencies (see Fig. 3,b) [33].
Finally, in Figure 11, the residuals of both derived GMMs (for PGA and PGV) are plotted with respect to magnitude, hypocentral distance and earthquake depth.The residuals do not follow any trend as a function of magnitude, distance or depth, which shows that the derived GMMs for the Gardanne data properly account for these predictor variables.Fig. 9. GMM for PGV (red lines) with two magnitudes (M w 0.5 and M w 1.0), compared with the data points used in the regression and with several GMMs from the literature.The dotted lines represent the þ/À 1 standard deviation boundaries.

Site amplification and site terms
Thanks to the adopted GMM formulation, it is also possible to estimate the station-to-station ground-motion site terms, for each of the 9 stations used in the dataset (see Table 3).The site term represents site-to-site variability [34] and it is useful to identify stations that tend to systematically under-or over-estimate the GMM predictions at these sites.These corrective terms may be seen as indicators of the sites associated with site amplification effects (which are not explained by the GMM).
According to a national level map of EC8 soil classes [35], the Gardanne area is mostly composed of EC8 soil class A (i.e., rock conditions).As a result, the locations of the recording sites are also assumed to correspond to rock conditions so that the derived GMM is considered to be mostly applicable to rock or stiff soil sites.The main consequence of this assumption is that the functional form of the proposed GMM (see Equation ( 4)) does not include any term related to site conditions, such as V s,30 -dependent factors.
However, a further study of potential site effects is also carried out in the Gardanne area, specifically at the location of the seismic stations.Horizontal to Vertical (H/V) spectrum ratios of recorded microtremors [36] at different stations of the seismic network are compared.For most stations, H/V ratios do not evidence any notable amplification effect (see example in Fig. 12 left), as expected.However, this is not the case for all stations, as some of them (e.g., SAVA, see Fig. 12 right) may be associated with significant site effects in the high-frequency range (15e20 Hz).
Also, some sites that were expected to be without any frequency amplification, as they are located on rock conditions, still generate high-frequency H/V peaks.These peaks are directional, as shown in Figure 13, but they are not created by a stable anthropic source.These high-frequency site effects, combined with the high-frequency excitation from small magnitude events, may explain why some small induced earthquakes may still be felt: as the small earthquake corner frequency is higher, there can be a non-negligible amount of amplified motion, thus reaching the perception threshold.
The estimation of statistical site terms through the GMM leads to the quantification of amplification effects at the recording stations.The positive site terms that have been found for SAVA station at short periods (see Table 3) seem to match the site effects in the high frequency range revealed by the H/V ratio.A similar observation may be formulated for the VERW station.On the other hand, no significant H/V peak has been found for the VILO station (see Fig. 12 left), and this is also confirmed by its site terms (mostly negative).

Shake-maps for recorded events
When an earthquake is detected (along with its characteristics, such as magnitude, epicentre location and focal depth), the application of a GMM provides estimates of the distribution of groundmotion parameters around the epicentre.If observations or records of the ground shaking (i.e., Fig. 10.Left: comparison of the response spectrum predicted by the GMM with a set of recorded spectra for a M w 1.0 event at a distance of 1 km (dashed lines represent þ/À 1 standard deviation); middle: comparison of the response spectrum predicted by the GMM with a set of recorded spectra for a M w 0.5 event at a distance of 1.5 km (dashed lines represent þ/À 1 standard deviation); right: normalized response spectra at T ¼ 0.1 s for increasing magnitudes, compared to a theoretical model of slope 2.
ground-motion measurements and macroseismic intensities) are available, they may be collected and used to update the distribution of the ground-motion field via statistical interpolation techniques.This process results in a shake-map [37], which provides ground-motion estimates constrained by field evidence expressed as intensity measures (IMs).
The review by Gu erin-Marthe et al. [38] of current shake-map algorithms and systems describes the ShakeMap® algorithm developed and operated by the U.S. Geological Survey [39], as well as alternative interpolation methods.It is proposed here to apply the Bayesian shake-map approach introduced by Gehl et al. [40], which is based on the updating of spatially correlated Gaussian fields.It is based on the following steps: Given the characteristics of the earthquake event, the prior distribution of the ground-motion field is provided by the application of the GMM.For each grid point i on the mapping area, thanks to the lognormal assumption used in most GMMs, a lognormal-normal conversion is able to express the conditional probability of Y i ¼ log PGA (or PGV) as a normal distribution, with the mean expressed as: where X i is the mean estimate of the logarithm of the ground-motion parameter from the GMM, f SS Fig. 12. Left: H/V spectral ratio at VILO station e rock soil; right: H/V spectral ratio at SAVA station.
Fig. 13.Left: H/V spectral ratio at VERW station; right: H/V spectral ratio at VERW station, with respect to azimuth.In the left plot, the norm of the total horizontal movement is used for each time window, while in the right plot, the direction-specific components of horizontal movement are averaged over all time windows.
is the standard deviation of the intra-event term, and t is the standard deviation of the inter-event term.The spatial correlation between intra-event terms is based on site-to-site distances and is estimated by the Jayaram & Baker [41]  The Bayesian method for the derivation of shakemaps is demonstrated here by applying it to one of the strongest events that occurred recently in the Gardanne area: April 19th 2019 event, M w 1.7 (estimated by Ineris), epicentre location at 5.5322 E Longitude and 43.4391 N Latitude, depth of 580 m (estimated by BRGM).
Ground-motion parameters (PGA, PGV and SA) have been recorded by the 9 monitoring stations for both events, as detailed in Table 4.
Once the data is collected, the following elements are used as inputs in order to generate the shakemaps: Characteristics of the earthquake event (magnitude M w , focal depth and epicentre location).Ground-motion parameters (PGA, PGV or SA) recorded by the monitoring stations.The Gardanne GMM for PGA, PGV or SA.The site terms (see Table 3) of the recording stations: these terms are subtracted from the logarithm of the recorded ground-motion parameters (see Table 4) in order to establish a baseline ground-motion field, without any site amplification factors.A map of site effects or amplification factors across the area: in the absence of any thorough study at the local level, a map of amplification established at the national scale [35] is applied here.As a result, the whole area around the Gardanne area is found to be associated with EC8 soil class A, with a soil amplification factor of 1 (i.e., no change from rock conditions).
The Bayesian shake-map delivers the outcomes within a couple of minutes, and the PGV and PGA shake-maps for the illustrative event are displayed in Figure 14.The contours of the estimated groundmotion fields are consistent with some of the values recorded by the nearby stations, e.g.PGA around 55 mg for the BULL station closest to the epicentre.Moreover, for this earthquake, two testimonies of macroseismic intensity III (macroseismic scale EMS-98 [42]) have been reported by citizens in the municipality of Gr easque.While the exact locations of these testimonies are unknown, the extent of Gr easque (south-east of the epicentre) is within the PGV range of 0.1e0.3cm/s, which is consistent with macroseismic intensity degree III, when applying, for instance, the ground-motion intensity conversion equation by Caprio et al. [43].
For this specific event, the residuals of GMM and shake-map estimates vs recorded values at the 9 stations are shown in Figure 15.The residuals of the initial GMM estimates are within the expected range of variability (i.e., around 0.2e0.3).When applying the specific site terms from Table 3, the range of the residuals becomes slightly more centered on the zero line, although this effect is not systematic for all stations.Then, after computing the shake-maps, all final residuals converge towards zero as expected: this result demonstrates the accuracy of the shake-map approach when constraining the predicted ground-motion field with respect to the observations.

Conclusions
This study has focused on the development of a GMM for post-mining induced seismicity, specifically based on ground-motion data collected from earthquakes in the Gardanne coal basin.One of the challenges has been the harmonization of a catalogue of recorded events, through the back-calculation of moment magnitudes with SourceSpec software.Potential differences between two catalogues (the operational one provided by Ineris and the one manually derived a posteriori by BRGM) and their respective estimates of the magnitude have been taken into account by assuming a standard deviation of 0.1 on the magnitude and by applying the Bayesian-based regression method by Kuehn & Abrahamson [14] for the derivation of the GMM.The formulation of the GMM has also led to the quantification of station-to-station ground-motion site terms: while the GMM is assumed to be derived for stiff soil or rock conditions, it is observed that some stations present positive site terms, which are confirmed by H/V ratios at high frequency on ambient noise measurements.When comparing the derived GMM with existing models for other types of induced seismicity (geothermal activity, hydraulic fracturing, and mining), it appears that the GMM by Cremen et al. (2020) provides results closest to the Gardanne specific model, although with a larger decay rate.A quadratic magnitude term has also been included in the model in order to better constrain predicted values at larger magnitude, in the case the GMM is used for larger earthquakes than the ones in the dataset (i.e., M w up to 1.7).Further developments could benefit from the collection of additional  datasets thanks to the seismic monitoring of other coal mine closures in Europe (e.g., [44]).Alternatively, analytical solutions based on a point-source model might also be combined with empirical data in order to derive a more robust model [45].
Finally, the derived GMM has been tested on a shake-map procedure in order to generate distributions of ground-motion fields in the Gardanne area, in case of an earthquake.It could, therefore, be integrated into the information bulletins released by Ineris.Such shake-maps prove useful to locate rapidly areas where potential damage may have occurred, as well as to obtain a reference map of the extent of the event, which may support further studies on the impact of post-mining induced seismicity on the population.

Fig. 1 .
Fig. 1.Situation map of the Gardanne coal basin (Fuveau-Gr easque area) and of the surrounding municipalities.

Fig. 2 .
Fig. 2. Left: Permanent seismic stations operated by Ineris (green triangles) and temporary seismic stations operated by BRGM (blue triangles).The black rectangle corresponds to the zoom area on the right; right: Epicenter location of 20 earthquakes (circles) between January 2019 and July 2019 analyzed in this study.The accelerometric network of BRGM (blue triangles) is used.
duration of 0.8 s, as shown in some examples in the Appendix (Fig.A1,b).FigureA1demonstrates an analysis of an earthquake (2019/05/02 09:13 M l ðBRGMÞ ¼ 0.36, M D ðInerisÞ ¼ 1.3, M w ðInerisÞ ¼ 0.1).Most of earthquakes occur within an epicentral distance of 2 km (focal depth less than 1 km), and the duration of the signals is very short.Therefore, the adopted time window (0.8 s) is sufficient to detect the characteristics of the waves.

Fig. 3 .
Fig. 3.The obtained source parameters for the 20 earthquakes (January 2019 to July 2019).(a) Moment magnitude (M w ) and local magnitude (M l ) in SourceSpec program.(b) Corner frequency (f c ). (c) Attenuation parameter (t * ).(d) Averaged hypocentral distance used for each earthquake.(e) Estimated Brune's stress drop (Dt) [21].(f) Estimated source radius.The error bars are also illustrated except for panel (d), in which the minimum and maximum values are indicated.

Fig. 4 .
Fig. 4. Comparison of magnitudes of the 20 earthquakes with respect to the obtained moment magnitude M w (Source_Spec).(a) Local magnitude obtained in SourceSpec, (b) local magnitude in BRGM catalog, (c) duration magnitude in Ineris Catalog and (d) moment magnitude in Ineris catalog.The solid and open marks present M w ðInerisÞ !0:3 and M w ðInerisÞ < 0:3, respectively.The broken lines show the best fitting line for the points of M w ðInerisÞ !0:3 with y ¼ Ax þ B, where x and y are the horizontal and vertical quantities of each panel.(a) A ¼ 0.770, B ¼ 0.987, (b) A ¼ 0.875, B ¼ 0.480, (c) A ¼ 1.065, B ¼ À0.689,(d) A ¼ 1.249, B ¼ À0.159.

Fig. 6 .
Fig.6.Map of the selected seismic events and recording stations for the construction of the GMM.Projection system EPSG:27573.

Fig. 11 .
Fig.11.Top: PGA residuals (GMM prediction PGA pred vs observed data PGA obs ) as a function of magnitude (left), hypocentral distance (middle) and earthquake depth (right); bottom: PGV residuals (GMM prediction PGV pred vs observed data PGV obs ) as a function of magnitude (left), hypocentral distance (middle) and earthquake depth (right).The dotted lines represent the þ/À 1 standard deviation boundaries.

Fig. 15 .
Fig.15.Residuals estimated for the event of April 19th 2019, regarding PGA (left) and PGV (right).GMM estimates are compared to recorded values without and with the site terms of the 9 stations (empty and full blue circles, respectively).Red circles represent the residuals of the shake-map estimates.

Table 1 .
Description of the stations currently operated by BRGM.

Table 2 .
Estimated coefficients for the GMMs in terms of PGA, PGV and SA.

Table 3 .
Estimated site terms for the 9 recording stations, in terms of PGA, PGV and SA.

Table 4 .
Recorded ground-motion parameters (PGA and PGV) for the April 19th 2019 event.