1 Introduction

The moment-based estimate of the stress drop is a source parameter that characterizes the amount of high-frequency energy radiated by an earthquake. For a given seismic moment, it is determined by the position of the corner frequency value and determines the amplitude of the high-frequency plateau in the acceleration Fourier amplitude spectrum (Brune 1970, 1971; Hanks and Johnson 1976). In recent decades, the stress drop computed for the Brune source model has attracted considerable interest due to its implications for various aspects of ground motion predictions (e.g., Baltay et al. 2013; Bindi et al. 2023a), the physical disparities between small and large earthquakes (e.g., Baltay et al. 2011), and frictional models of earthquake rupture (e.g., Abercrombie 1995). The stress drop is also related to one of the open questions in seismology, namely whether the scaling of the stress drop is magnitude-dependent or self-similar. For moderate to large earthquakes, the magnitude-independent or self-similar scaling of the stress drop is supported by theoretical (Aki 1967; Haskell 1966) and observational considerations (Kanamori and Anderson 1975). However, the breakdown of self-similar scaling has been observed at both regional and global scales, especially for smaller-magnitude events (e.g., Bora et al. 2013; Bindi et al. 2020; Kanamori et al. 1993; Oth et al. 2011; Trugman et al. 2017).

The discrepancy between stress drop estimates obtained by different methods may further increase the controversy (e.g., Abercrombie 2014; Kaneko and Shearer 2015; Neely et al. 2020; Shearer et al. 2019; Kemna et al. 2021). Different methods for estimating stress drop have their own assumptions and limitations (e.g., different data types, and data processing) that can lead to conflicting results and interpretations. For example, estimating stress drop from the corner frequency requires many assumptions to be made in setting up the source model, such as the shape of the rupture area and the average rupture velocity, which makes it difficult to compare different studies as these assumptions often vary. In addition, a wide frequency bandwidth is required in the recorded data to estimate corner frequencies over a large range of magnitudes.

In the last few decades, the expansion of seismic networks at both regional and global scales, as well as the establishment of a standard format for archiving and exchanging waveforms and metadata by the International Federation of Digital Seismograph Networks (FDSN), has led to an exponential increase in available data. Recent developments in processing software (e.g., Zaccarelli et al. 2019) have also enabled researchers to efficiently process large amounts of waveform data. We took advantage of these developments to achieve our goal of systematically estimating the stress drop from the harmonized dataset at the European scale. For this purpose, we downloaded a large set of available records of seismic events in Europe (mainly in Germany, Italy, France, and Poland) with magnitudes between 1.5 and 6.5 at distances less than 200 km through the European Integrated Data Archive, EIDA (Strollo et al. 2021), and developed a procedure to detect outliers and extract high-quality records.

In this study, we isolate the source contribution to ground shaking from propagation and site effects by performing a spectral decomposition (Oth et al. 2011). To account for significant attenuation differences within the study area, we first calibrate the spectral attenuation models for two separate regions, mainly dividing the study area along the Alps, which is consistent to the regionalization introduced by previous studies (e.g., Bindi et al. 2019; Bindi and Kotha 2020; Weatherill et al. 2020). The obtained source spectra are fitted with a standard ω2-model (Brune 1970, 1971) to estimate the seismic moment and the corner frequency and, in turn, to compute the stress drop. Finally, we discuss the scaling of the stress drop with the seismic moment obtained for the whole study under uniform processing and model assumptions.

2 Data processing and selection

The data download process is guided by a seismic event catalog from the event web service of the International Seismological Centre (ISC, http://www.isc.ac.uk/fdsnws/event/1/query). We select events with focal depths shallower than 60 km to include the crustal events and with magnitudes between 1 and 6.5 between January 1990 and May 2020. The total number of the available three-component records in our database is about 18 million waveforms recorded by high sampling rate channels (i.e., velocimetric channels HH and EH, and accelerometric channels HN, HL, HG).

To ensure that events from different hypocentral distances are fully captured in the target windows, recordings are 4 min long starting 1 min before the theoretical P-wave arrival time (computed using the AK135 global velocity model). As the download window of the signal is large and it could contain multiple events and potential outliers, multiple event detection, outlier detection, and quality control algorithms are applied in the data processing. The stream2segment software (Zaccarelli et al. 2019) was used to process the recordings following previous studies (e.g., Bindi et al. 2019). The instrumental response is removed in the spectral domain, with the water level regularization set to 100 dB, and the processed recordings are converted to acceleration. The bandpass filter is designed as an a-causal Butterworth filter with a high-pass cutoff frequency of 0.5 for M < 3, 0.3 for 3 < M < 6, and 0.08 Hz for M > 6, respectively. The low-pass cutoff frequency is set at 40 Hz. The signal windows correspond to the interval between the 2.5 and 97.5% of the cumulative squared velocity computed from the estimated P-wave arrival, tapered to zero at both ends with a 5% cosine profile. The frequency range of interest is from 0.3 to 25 Hz with 48 frequency points.

In order to obtain a high-quality dataset for the spectral decomposition analysis, the signal-to-noise ratio (SNR) is calculated with respect to the pre-event noise over the bandwidth [0.3, 25] Hz, and only recording with SNR > 20 is retained for the analysis. In addition, outliers are identified and removed using the modified Z-score method. This approach is applied to the residuals computed with respect to the Ground Motion Prediction Equation (GMPE) of Bindi et al. (2014), considering peak ground velocity (PGV) and peak ground acceleration (PGA), with the threshold score for identifying outliers set to 3.5 in absolute value (Iglewicz and Hoaglin 1993). Finally, only earthquakes and stations with at least 3 recordings are considered for the spectral decomposition.

The final dataset contains 220,000 spectra, calculated as the square root of the sum of the two horizontal components squared. These recordings are associated with more than 6500 earthquakes recorded by a total of about 1600 station channels. The distribution in terms of geographical location, hypocentral distance, and magnitude of the events and stations included in the spectral analysis is shown in Fig. 1. To account for regional propagation patterns which are consistent with the models derived from previous studies (i.e., Bindi et al. 2019; Bindi and Kotha 2020), the data are grouped into two regions according to the location of the events, hereafter referred to as A (Central Europe and Balkans) and B (mainly Italy). For the final source spectrum fitting, 6227 spectra with a hypocentral distance of less than 150 km are selected due to the stability of the resulting site and source terms, which will be discussed in the following section. The selected spectra contain 98% of natural tectonic events and 2% of non-tectonic events. The events of the selected spectra are matched against the induced seismicity catalog (Grünthal 2014) to identify the non-tectonic events (potentially induced events) in our dataset, which include geothermal projects, hydrocarbon exploration, and resource mining.

Fig. 1
figure 1

a Map showing earthquakes (circles) and stations (open squares) used for the spectral decomposition for regions A and B. Circle colors indicate the regionalization. The black box indicates the transition zone we selected for the further attenuation check. b Catalog magnitude obtained from ISC through the FDSN event query (mostly mb and ML) versus the number of segments. c Hypocentral distance versus the number of segments. The black line indicates the maximum distance considered for fitting source spectra

Full size image

3 Methodology

3.1 Spectral decomposition

To isolate the source, propagation, and site contributions to the ground motions in the frequency domain, we perform a non-parametric spectral decomposition approach known as Generalized Inversion Technique (GIT; Castro et al. 1990). The Fourier Amplitude Spectra (FAS) at frequency \(f\) of the recordings are represented as logarithms of the source spectrum \({S}_{i}\left(f\right)\), the attenuation along the travel path \(A\left(f,r\right)\), where r is the hypocentral distance, and the site amplification \({Z}_{j}\left(f\right)\). Considering the sequences of earthquakes with sources \({S}_{i}\left(f\right)\), with \(i\) = 1, …, N, recorded by stations characterized by the site responses \({Z}_{j}\left(f\right)\), with \(j\) = 1, …, M, and the regionalized attenuation models \({A}^{k}\left(f,{r}_{ij}\right)\), the equation can be written as

$${{\text{log}}}_{10}{U}_{ij}\left(f,{r}_{ij}\right)={{\text{log}}}_{10}{S}_{i}\left(f\right)+{{\text{log}}}_{10}{A}^{k}\left(f,{r}_{ij}\right)+{{\text{log}}}_{10}{Z}_{j}\left(f\right)$$
(1)

with k = A, B identifying two different spatial domains for the attenuation models.

To assess the stability of the inversion results, we perform 100 bootstrap inversions at each frequency point following the procedure described in Parolai et al. (2000). A smooth function of distance, \(A\left(f,{r}_{ij}\right)\) is constrained to be one at 10 km. The distance is subdivided into distance bins of 5 km width for the determination of the attenuation functions. To capture the regional patterns of propagation in a large study area, the recordings are divided into regions A and B (Fig. 1).

In this study, the decomposition in Eq. (1) is solved twice:

  1. (1)

    First, we focus on the regional attenuation models by solving Eq. (1) in a least-squares sense, considering both source, propagation, and site terms (single-step decomposition). Since the resulting attenuation functions from the single-step decomposition are independent of the reference site condition (Oth et al. 2011), we constrain the average site response of all stations to one regardless of the frequency. It is noted that to derive the stable attenuation models, the hypocental distance is set up to 200 km.

  2. (2)

    The second decomposition is performed considering the FAS corrected for the attenuation term provided by the first decomposition (Eq. (2)). The aim of the second decomposition is to separate site and source contributions by a priori constraining the site response of a selected reference station. Since our target is to characterize the source parameters, we tried to cut the distance as short as possible to keep only the source effects. We select the subsets with different distances (i.e., 60, 100, and 150 km) to test the trade-off between the concentration of the source effect and the overlapping path between stations and events (Figure S1). The resulting site amplifications with distance < 150 km are a lower site-to-site variability of the site response and source spectra shape which are more ω2-shaped. Regarding the trade-off between the concentration of the source effect and the overlapping path between stations and events, we consider the hypocentral distance of 150 km.

To prepare the input dataset for the second part of the decomposition, the FAS are corrected for the regionalized attenuation models which have been derived in the first decomposition as

$${{\text{log}}}_{10}{R}_{ij}\left(f\right)={{{\text{log}}}_{10}{U}_{ij}\left(f,{r}_{ij}\right)-{\text{log}}}_{10}{A}^{k}\left(f,{r}_{ij}\right)$$
(2)

and the FAS residuals \({{\text{log}}}_{10}{R}_{ij}\left(f\right)\) are used as input for the second part of the decomposition as

$${{\text{log}}}_{10}{R}_{ij}\left(f\right)={{\text{log}}}_{10}{S}_{i}\left(f\right)+{{\text{log}}}_{10}{Z}_{j}\left(f\right)$$
(3)

The second part of the decomposition (Eq. (3)) allows the separation of source and site effects. To determine site response functions and source spectra relative to the same reference condition, and to remove the trade-off between the source and site terms, we defined a reference site at station CH.LLS (station Linth-Limmern), which is installed on a rock with shear-wave velocity averaged over the uppermost 30 m (VS30) of 2925 km/s (Fäh et al. 2009). The site amplification at the reference station is constrained to a function equal to 1 for frequencies below 10 Hz and a function of \({e}^{-\pi {\kappa }_{0}\left(f-10\right)}\), with \({\kappa }_{0}\)= 0.007 s (Pilz et al. 2019) above 10 Hz, to account for near-surface attenuation effects at high frequencies (Anderson and Hough 1984).

3.2 Spectral fitting and determination of stress drop

The spectral fitting was carried out on the inverted source spectra to estimate corner frequency \({f}_{C}\) and seismic moment \({M}_{0}\) following Eq. (4). The latter consists of a standard ω2-model multiplied by an exponential factor, \({\kappa }_{source}\), applied to frequencies above \({f}_{k}\):

$$\begin{array}{cc}S\left(f\right)={\left(2\pi f\right)}^{2}\frac{V{R}^{\theta \varnothing }F}{4\pi \rho {{v}_{S}}^{3}{R}_{0}}\dot{M}\left(f\right){e}^{-\pi {\kappa }_{source}\left(f-{f}_{k}\right)}, with\;\dot{M}\left(f\right)=\frac{{M}_{0}}{1+{\left(\frac{f}{{f}_{C}}\right)}^{2}}& ,\end{array}$$
(4)

where \(S\left(f\right)\) represents the acceleration source spectrum at the reference distance \({R}_{0}\) = 10 km, \(\dot{M}\left(f\right)\) denotes the moment-rate spectrum, \(V=\frac{1}{\sqrt{2}}\) is a factor to account for the partition of total shear-wave energy into two horizontal components, \({R}^{\theta \varnothing }\) is the average radiation pattern of S-waves set to 0.55 (Boore and Boatwright 1984), \(F\) = 2 is the free surface factor, \(\rho\) = 2.8 g/cm3 is the density, and \({v}_{S}\) = 3.3 km/s is the shear-wave velocity near the source.

Stress drop and seismic moment are two key parameters of stochastic, hybrid, and physics-based earthquake rupture and ground motion models (e.g., Razafindrakoto et al. 2021). To ensure that the earthquake magnitudes and stress drops resulting from our study are consistent with the seismicity and ground motion models of past, present, and future European seismic hazard studies. We first perform a preliminary fit to estimate the average bias of the seismic moments with respect to the moment magnitude of the European Mediterranean Earthquake Catalogue (EMEC; Grünthal and Wahlström 2012) which is the moment magnitude scale used in European seismic hazard studies (Danciu et al. 2021). The average bias is computed considering earthquakes in the range of M3.5–6.3. All source spectra are shifted by an average of 0.5 magnitude units. The final corrected spectra are re-fitted again to a ω2-model and the magnitude is constrained to the EMEC magnitude for the M > 5 events, if the corresponding magnitude is available in the EMEC.

The obtained corner frequency and moment magnitude are used to compute the stress drop \(\Delta \sigma\), assuming a circular rupture model with uniform stress drop as Eq. (5) (Hanks and Thatcher 1972):

$$\Delta \sigma =8.5{M}_{0}{\left(\frac{{f}_{C}}{{v}_{S}}\right)}^{3}$$
(5)

where \({v}_{S}\) = 3.3 km/s for all events (Lu et al. 2018).

4 Results and discussion

4.1 Regionalized attenuation models

The attenuation models for regions A (Fig. 2c) and B (Fig. 2d) show similar patterns below 30 km whereas, while at distances above about 50 km, the attenuation becomes stronger in region B, especially at frequencies > 10 Hz, in agreement with previous studies (e.g., Bindi and Kotha 2020). Both attenuation models indicate a deceleration in attenuation between 70 and 100 km, attributed to phenomena such as Moho reflections and intra-crustal discontinuities.

Fig. 2
figure 2

Attenuation models versus hypocentral distance for a the no regionalization (the whole study area), b the transition (the area between regions A and B), c region A, and d region B. The colored curves represent the attenuation for different frequency ranges. Black dash lines indicate attenuation rates proportional to the inverse of the reference distance and the inverse of the square of the reference distance

Full size image

To check the effect of regionalization on the attenuation models, we repeat the decomposition but either introduce no regionalization (i.e., regions A and B are merged) or consider a small transition zone that partially overlaps regions A and B (black box in Fig. 1). Since the non-parametric decomposition is a data-driven approach, the attenuation for the case of no regionalization (Fig. 2a) is similar to the attenuation obtained for region B, where most of the earthquakes are located. The transition zone contains 9083 records from region A and 8361 records from region B. As the dataset for the transition zone is well balanced in terms of the contributions from regions A and B, the obtained spectral attenuation (Fig. 2b) captures the transition from weaker to stronger attenuation as one moves from region A to B. In the following, we consider the attenuation models for regions A and B to correct the FAS for large-scale attenuation effects.

4.2 Non-parametric site amplifications and source spectra

The resulting site amplification at all considered stations is the amplification relative to that of the reference site, CH.LLS (Fig. 3 and the electronic supplement). However, the amplification at the NRCA station (Norcia, IV network) is well captured in this study, which has a peak at around 7 Hz as found in other studies as well (Vassallo et al. 2022). It is noted that the number of recordings at some stations is too low (n = 3) to catch the stable site response, and also the seismic instrument filter affects the amplification above 20 Hz. The stability and reliability of the decomposition results still have to be evaluated since our study area is relatively large and contains various site conditions.

Fig. 3
figure 3

Site amplifications at all considered stations. The red solid curve represents the site amplification of the reference site, CH.LLS. The blue dashed curve indicates the site amplification of the station NRCA in the IV network

Full size image

Therefore, we compare site amplifications with those derived in Bindi et al. (2023a) which are constrained to the same reference site. A selection of site amplification at stations is shown in Fig. 4. The amplitude at station CH.LLS at low frequencies (< 1 Hz) and some stations (FR.OGS1 and IV.NRCA) at higher frequencies (20 Hz) has a factor of 2 difference between different channels, which might be due to the difference in the data selection and the data processing. However, the patterns and amplification values are in good agreement with those of Bindi et al. (2023a) over most stations.

Fig. 4
figure 4

Comparison of the ratio of site amplifications to the reference site. Solid curves indicate the site amplification from this study. Dashed curves indicate the site amplification from Bindi et al. (2023a)

Full size image

The obtained source spectra with regionalized attenuation corrections are shown in Fig. 5. Since the source spectra are derived in the same spectral decomposition with the same reference site, it is possible to directly compare the results of large events. As expected, the spectral amplitudes of significant events are higher for the larger events, such as the 2009 L’ Aquila, 2012 Emilia, and 2016 Amatrice, Norcia, Visso earthquakes.

Fig. 5
figure 5

Acceleration source spectra derived from the GIT inversion (entire study area). The source spectra with the regionalized attenuation correction for a region A (north) and b region B (south). Colored curves show the selected events (the 2009 L’ Aquila, 2012 Emilia, and 2016 Amatrice, Norcia, Visso earthquakes)

Full size image

4.3 Regionalization impact on source spectra and stress drop

In order to assess the impact of the regionalized attenuation models on the source spectra, we compare the ratio of spectra between no-regionalized attenuation and regionalized attenuation cases (Fig. 6). For region A, the amplitudes of the source spectra at 10–25 Hz with regionalized attenuation corrections are generally lower than those with no-regionalized attenuation corrections. The maximum difference is a factor of 4. For region B, the amplitudes with regionalized attenuation corrections are higher but the difference between the two cases is smaller, within a factor of 2. Therefore, the regionalization systematically affects at the high frequencies in both regions. It is again emphasized that the regionalization is essential to account for first-order attenuation variations over such a large region of interest.

Fig. 6
figure 6

The ratio between spectra derived from the no-regionalized attenuation correction and the regionalized attenuation correction for a region A and b region B. Colored curves show the selected events that occurred in Italy (the 2009 L’ Aquila, 2012 Emilia, and 2016 Amatrice, Norcia, Visso earthquakes)

Full size image

Once the corrected source spectra are obtained, the source spectra are fit with a ω2-model to derive the corner frequency of each spectrum. The corner frequencies and moment magnitudes obtained are shown in Fig. 7. Due to the limited observational bandwidth limitation of surface stations, the high-frequency energy (20–25 Hz) is depleted, biasing the measured and estimated \({f}_{c}\) to lower values, strongly affecting smaller earthquakes with higher theoretically expected \({f}_{c}\) values, as shown in Kemna et al. (2021). Consequently, we observed that the obtained corner frequency did not increase for decreasing magnitudes (M < 3). Thus, the robust resulting magnitude and the corner frequency range in this study are M > 3 and 0.3 to 20 Hz. According to the standard deviations of the fitted corner frequencies and moment magnitudes, most of the source spectra follow the standard ω2-model, and their corresponding corner frequencies and moment magnitudes are reasonable in the effective area (the gray area in Fig. 7). However, about 1% of the total spectra do not fit with the ω2-model or have a high standard deviation. Fifty-seven percent of these spectra are non-tectonic events in Groningen and Poland and 43% are small earthquakes in southeastern Italy. We attribute this discrepancy to the limitations of our large-scale dataset to perform spectral decomposition in regions that are not sampled densely enough. The limited redundancy in the dataset for some areas, such as those in central and northern Europe, does not allow the source term to be properly isolated from site and attenuation contributions. This is the case, for example, for the non-tectonic events occurring in Groningen, for which we obtained unrealistic source spectra, whereas a recent spectral decomposition study performed with a local dataset (Ameri et al. 2020) obtained source spectra well described by a Brune ω2-model. This suggests the need for focused analysis and data collection for these small and specific events.

Fig. 7
figure 7

Scaling of corner frequency versus moment magnitude obtained by fitting the source spectra with a Brune ω2-model. Blue circles with bars indicate the results obtained with the regionalized attenuation correction. Orange circles with bars indicate the results obtained with the no-regionalized attenuation correction. Gray circles indicate the events for which the standard deviations are large or for which the spectra do not fit with the ω2-model

Full size image

When fitting the source spectra and deriving stress drops, we also observed a strong impact of regionalization. On the one hand, the obtained corner frequencies show that the regionalized attenuation model makes the corner frequency values more variable (Fig. 7). On the other hand, the ratio of the resulting stress drops between the regionalized and the no-regionalized attenuation model shows that the stress drops are lower in region A and higher in region B when the regionalized attenuation model is considered (Fig. 8a). Accordingly, if we correct the spectra with the no-regionalized attenuation model, the stress drops obtained may be overestimated in region A and underestimated in region B. It is then obvious that considering regional variations in the attenuation affects the quantification of stress drops. Therefore, in the following discussion, we only consider stress drop values resulting from regionally dependent attenuation models. Figure 8b shows the geographical distribution of the stress drops. Stress drops in region A are generally smaller than those in region B. However, it should be noted that the event magnitudes in region A are mostly smaller than those in region B. The geographical distribution, which mixes all earthquake magnitudes, would bias the interpretation of the regional variability and magnitude dependence of stress drops. Therefore, in the next section, we divide the events into magnitude bins to discuss the stress drop probability distributions fairly.

Fig. 8
figure 8

a Map of the ratio of the stress drops between without and with regionalized attenuation models. b Map of the stress drop with regionalized attenuation models

Full size image

4.4 Stress drop probability distributions

As we mentioned in Sect. 1, whether the stress drop is magnitude-dependent or self-similar scaling is a matter of debate, and it is difficult to compare between different studies due to a number of assumptions made in each study. In this study, we are able to scale the stress drops for a large number of small to moderate earthquakes under the same assumptions. Figure 9a shows the scaling relationship of stress drops with the moment magnitude for all events. The moment magnitudes are mainly between 2.5 and 3.5, and the mean stress drop of the magnitude larger than 3 is 8.63 MPa with a log standard deviation of 1.59. We find a positive correlation between stress drop and earthquake magnitude for small magnitudes ranging between 3 and 4, and a self-similarity for events of magnitude larger than magnitude 4 (mean stress drop of 18.25 MPa), which is consistent with previous studies. Indeed, previous studies have found a constant stress drop above magnitude 4 in southern Europe, the Mediterranean region, and the Middle East (Bora et al. 2017; Edwards and Fäh 2013). In addition, Allmann and Shearer (2009) showed that the median stress drop of about 4 MPa is independent of the seismic moment, implying self-similarity over the moment magnitude ranging 5.2 to 8.3 range in their global dataset.

Fig. 9
figure 9

Scaling relationship of a stress drops versus moment magnitude (all events in the study area) and b stress drops versus depth. Red vertical bars are the mean ± one standard deviation. Orange diamonds indicate the stress drops derived by Bora et al. (2017). Blue diamonds indicate the stress drops derived by Edwards and Fäh (2013). The dashed line in a indicates the stress drops associated with a corner frequency of 25 Hz. c Stress drops versus moment magnitudes for the two regions considered in this study. d Stress drops versus moment magnitude according to the source origin (tectonic or non-tectonic)

Full size image

Whether stress drop is depth-dependent in the crust remains a question in the community due to the large uncertainties in stress drop measurements; some studies have found that stress drop increases with depth (e.g., Boyd et al. 2017; Huang et al. 2017) and others have not (e.g., Allmann and Shearer 2009; Shearer et al. 2006). Recently, Abercrombie et al. (2021) concluded that the depth dependence of the stress drop systematically decreases when a depth-dependent attenuation model is considered. On the other hand, Bindi et al. (2023b) showed that the tendency of the stress drops to increase with depth is reduced when depth-dependent attenuation and velocities are considered. In our study, we observed that the average stress drop of each magnitude bin increased with depth (Fig. 9b), similar to other studies where this trend was observed when no depth-dependent attenuation model was considered. However, the lack of accurate focal depth and the few deeper events would bias the observed trend of stress drop with a depth, and we therefore refrain from further interpreting this depth dependence in physical terms.

Although some of the non-tectonic events do not fit to the ω2-model well, the remaining ω2-shape consistent non-tectonic events still provide an opportunity to shed some light on the variation of the stress drop with respect to the source origin at first glance. We observe that the stress drop of non-tectonic events appears to be self-similar instead of following the scaling relationship with the moment magnitude of the tectonic events (Fig. 9d). However, the limited number of non-tectonic events in our database and the assumptions for their source model make the uncertainties very large; therefore, we suggest including more non-tectonic events in future analysis to clarify the scaling relationship for the non-tectonic events.

Many methods for simulating ground motion from earthquakes have their particular strengths, but they depend, directly or indirectly, on the selected stress drop and its distribution as input parameters. The distributions of stress drop values and their regional dependencies are, therefore, important for the calibration of these simulations. The stress drops found in this study show significant regional variations (Fig. 9c). The stress drops are mainly distributed in the range of 0.01–10 MPa in region A (a mean of 0.62 MPa) with a magnitude range of M2–3 and 0.1–100 MPa in region B (a mean of 8.08 MPa) with a magnitude range of M2.5–3.5. Since there is an uneven number of earthquakes in different magnitude bins, we further split them into a subset of M < 3, 3 ≤ M < 3.5, 3.5 ≤ M < 4, and M ≥ 4 for both regions (Fig. 10). The stress drop probability distributions in different magnitude bins show that the stress drops in region B are generally higher than those in region A for the same magnitude bin. The stress drops in region B (south) increase with earthquake size, whereas the stress drops in region A (north) are more likely to be self-similar.

Fig. 10
figure 10

Stress drop probability distributions according to region and magnitude bins

Full size image

The mean value of the stress drop is a critical factor in calibrating the level of high-frequency shaking for ground motion simulations. Cotton et al. (2013) suggested that the variability of the stress drop as an input to ground motion simulations is an equally important factor, with the value having a large effect on the variability of the simulated ground motion. The stress drop variability from other spectral studies and our results is between 0.8 and 1.8 (Table 1). In our study, we also estimate the variability of stress drops for each magnitude bin in both regions and the entire study area listed in Table 2. The stress drop variability is lower in region B, but the variability is quite high in region A and in the whole region. These regional differences remain for each magnitude bin.

Table 1 Stress drops and variability from other spectral studies (\({M}_{0}{\text{,} \, {\text{f}}}_{c}\))
Full size table
Table 2 Stress drops mean value and standard deviation per magnitude bin and region
Full size table

5 Conclusions

The aim of this study is to investigate the characteristics of the source parameters and the scaling relationship of the stress drop with moment magnitude for earthquakes occurring in central and southern Europe. We apply a non-parametric spectral decomposition approach to a large dataset of recordings retrieved from the EIDA archive. We introduce a regionalization for the attenuation models considering two regions divided along the Alps, and show the significant impact of the regionalization on the resulting source spectra and stress drops.

To fit the standard ω2-model, we merged the records from both regions for the site and source decomposition, and constrained the site amplification to a single reference site, CH.LLS, to derive comparable source spectra. Most of the source spectra obtained follow the standard ω2-model and a corner frequency and moment magnitude can be estimated. About 1% of the total spectra cannot be fitted with a ω2-model: 57% of these are non-tectonic events in the Groningen region and Poland, and 43% are small earthquakes in southeastern Italy. This discrepancy is due to the limitations of our large dataset to perform spectral decomposition in regions that are not sampled densely enough. The limited redundancy in the dataset for some areas (e.g., Groningen) does not allow the source term to be properly isolated from site and attenuation contributions. This suggests the need for focused analysis and data collection for these small and specific events.

We find a positive correlation between stress drops and earthquake magnitudes in the magnitude range of 3 to 4, and a self-similarity for events with a magnitude larger than 4. The mean stress drop (13.8 MPa) is consistent with other studies using European and global data. A regional variation of stress drops was also found: stress drops are generally higher in region B (south, more active) than in region A (north). In the southern part, the stress drops increase with earthquake size, whereas, in the northern part, the stress drops are more likely to be self-similar.

The results of this non-parametric spectral decomposition demonstrate the potential of working on a large dataset and considering the regionalized attenuation model, and also provide probability distributions for regional European earthquake stress drops that can be used to calibrate stochastic, hybrid, or even physics-based simulations (e.g., Bindi et al. 2023a).