Statistical Analysis of Mt. Vesuvius Earthquakes Highlights Pitfalls in Magnitude Estimation

Abstract

Here, we characterize the statistical behaviour of the Mt. Vesuvius seismicity using distinct available catalogues. Our analysis confirms that for this area, the GR distribution exhibited two scaling regimes of the b-value, not commonly observed for the standard frequency-magnitude distribution of earthquakes. By assuming a physical cause, we tested four different hypotheses for the source of the break in the scaling: finite size effect, depth variations in the b-value, radial dependence in the b-value, and different b-values for swarm and non-swarm events. None of the above reasons are able to explain the observation. Thus, we investigated the possibility of some pitfalls in magnitude estimation. Based on our analysis, we suggest there is a bias in the duration magnitude the catalogues are based on. This is due to the arbitrary extrapolation to smaller magnitudes of a linear regression derived for earthquakes with $m\ge 3.0$. When a suitable correction is applied to the estimated magnitude, the GR distribution assumes the usual shape, with a b-value closer to that usually observed in volcanic areas. Finally, the analysis of the time variation of some statistical parameters reveals that the state of the volcano appears to be stationary over the entire analysed period, possibly with only a slight increase in the b-value, indicating a small reduction in differential stress.

1. Introduction

Mt. Vesuvius is a composite strato-volcano located in the Campania Plain (Southern Italy), belonging to the Somma–Vesuvius volcanic complex. It is 1281 m high and 10 km wide. Some eruptive fissures, aligned in the E–W and N–S directions, cross the Mt. Somma caldera and the southern and western flanks of Vesuvius [1]. A melted zone has been individuated at ∼8 km depth [2,3], whereas the presence of a deep reservoir at depths of about 2.5–5.0 km b.s.l. has been hypothesized by Chiodini et al. [4], Del Pezzo et al. [5].

The seismicity of Mt. Vesuvius is characterized by low energy earthquakes, mainly distributed around the crater axis, at depths between 0 and 1 km (Figure 1). The strongest event since the end of the last eruption in 1944 occurred on 9 October 1999, with a magnitude of $m_d$ = 3.6 [6]. The regional tectonic stress and the gravitational load—in association with both the hydrothermal activity in the crater area and episodes of injection of deeper magmatic fluids—should be at the basis of the local seismicity [4,7].

The b parameter of the Gutenberg–Richter (GR) distribution [8] represents one of the most important parameters in the characterization of seismic occurrence because (i) it is crucial in the evaluation of the seismic hazard and (ii) it is inversely correlated to the stress state [9,10,11,12].

The spatial distribution of the b-value has been largely used to characterize stress regimes [13] and the spatial variations of the stress intensity in different seismogenic areas [14,15,16], including intermediate and deep seismic zones [17]. Several studies have also investigated the distribution of the b-value to individuate asperities on active faults [14,15,18], to characterize the complexity of the fault geometry [19], or to discriminate between distributed and on-fault seismicity [20]. Its analysis can also help in enlightening the presence of magmatic chambers in volcanic areas [15,21,22,23].

The time variations in the b-value represent a useful instrument for characterizing the nature of earthquakes in seismic sequences. Indeed, the variation in the b-value with time has been investigated to highlight the differences between foreshocks and background seismic activity [24,25,26,27,28,29] or to predictively discriminate between large foreshocks and mainshock [30]. Also, by stacking many earthquake sequences, an increase in b-value has been observed right after the occurrence of a mainshock [31].

Differently from tectonic earthquakes, quite often, volcanic events cluster in space and time without a triggering event of a higher magnitude [32]. These kinds of seismic sequences are called swarms. In general, they are recorded in either volcanic or geothermal areas, or even in tectonic environments [33,34,35]. It has been suggested that swarms occur in regions characterized by a high heterogeneity in terms of their material properties and stress concentration [36]. Their activity has been associated with stress changes induced by aseismic processes such as pore pressure changes [37] or fluid intrusion [35,38,39,40,41]. Volcanic swarms are usually the main reported seismic precursor for volcanic eruptions, especially for volcanoes that have been silent for decades or longer [33]. Although many models have modified the original ETAS model [42,43] and been able to reproduce some statistical feature of the seismic swarms [44,45,46,47], they have not been able to fit the occurrence rate: neither the Omori law [48] nor a simple relationship can describe the temporal evolution of the volcanic earthquake’s swarms well. The reason for such a difficulty is the duration of the swarms, which are often very short and do not provide a sufficient number of events for a reliable statistical analysis. However, by stacking many swarms with an average rate of occurrence, an analytic expression that reproduces the experimental observations can be derived [49].

The inter-event time between two successive earthquakes, $\mathsf{\Delta }t$, represents a powerful instrument to investigate the clustering properties of seismic events [50]. The $\mathsf{\Delta }t$ distribution $p\left(\mathsf{\Delta }t\right)$ can be considered universal when the inter-event times $\mathsf{\Delta }t$ are rescaled by the mean occurrence rate, R [51,52,53,54]. In other words, R defines a ‘local’ time scale that characterizes the earthquake occurrence, whereas their clustering properties can be considered universal. This result was firstly obtained for pseudo-stationary periods, revealing that earthquakes tend to cluster even if their occurrence is apparently Poissonian [52]. The universality of $p\left(\mathsf{\Delta }t\right)$ has also been observed for non-stationary periods [55] and for aftershock sequences [56,57]. Although the universal behaviour of $p\left(\mathsf{\Delta }t\right)$ has been questioned [58,59,60,61,62,63,64], the departure from universality, observed at small $\mathsf{\Delta }t$, has been explained in terms of four typical time scales [65]: the inverse rate of independent events, $\lambda$; the mean inverse rate of correlated events; the time parameter c defined in the Omori law [48]; and the catalogue duration T (this latter is irrelevant for the analysis presented here). A simpler approach defines a new expression for $p\left(\mathsf{\Delta }t\right)$, including the departure from universality [66]. Exception to this behaviour is observed for volcanic explosion-quakes (i.e., Stromboli), which may exhibit a Poissonian behaviour [67] that can be explained in terms of bubble formation in the magma chamber.

Here, we analyse the seismicity of mount Vesuvius, evidencing some peculiar behaviour of the Gutenberg–Richter (GR) distribution.

2. Data: The Catalogues of the Mt. Vesuvius Instrumental Seismicity

Although the settlement of the Vesuvius Observatory (National Institute of Geophysics and Volcanology; INGV-OV) dates back to 1841 and the detection of local seismicity has been carried out since the first seismic detector was installed by [68], a regular collection of data relative to the occurrence of earthquakes has not been preserved up until about the half of the last century (e.g., [69,70,71]). However, because of discontinuity in both the number and type of the installed instruments, the data cannot be considered homogeneous, at least for some further decades.

At present, three catalogues of the Mt. Vesuvius instrumental seismicity are available for previous decades, starting in different years (

Table 1. Mt. Vesuvius instrumental seismicity catalogues. Only events with an assigned magnitude are considered.

Catalogue	Starting Time	Ending Time	No. of Earthquakes
OVO	23 February 1972	27 December 2021	11,679
BKE	1 January 1999	31 December 2021	19,096
ALL	1 January 1999	8 September 2023	20,615
LOC	1 January 1999	8 September 2023	8766

). Two of these were based on single-station detection and magnitude estimation (OVO, BKE; Figure 2), these databases [72,73] cover up to 2021 and, according to the authors, they will be updated periodically, supposedly for historical continuity. The third one (ALL) is the current ordinary catalogue. It relies on the recordings of the monitoring seismic network, run by INGV-OV, and is regularly compiled according to modern standards.

The OVO and BKE catalogues cover the periods 1972–2021 and 1999–2021 (Figure 3), respectively. For both databases, each record contains the date, origin time, and event magnitude. The latter is only available as the duration magnitude $m_d$, which is the only one usually estimated at INGV-Osservatorio Vesuviano. The $m_d$ at OVO is determined from the coda duration, $\tau$, and is evaluated manually by the operator and converted according to the formula:

$$m_{d_{OVO}}=2.75log\tau -2.35.$$

(1)

This equation was derived by analysing regional earthquakes and comparing the waveform duration observed at OVO with the local magnitude $m_L$ for the same event, as determined from the recording of the Wood–Anderson seismometer installed at Roma Monte Porzio (RMP), about 180 km apart [74]. Because of the distance between the two stations, the comparison was limited to events with $m\ge 3.0$, then extrapolated to lower magnitudes. Later on, starting from this parent relation, an analogous relation was obtained for BKE by comparing the duration observed at this station with the one determined at OVO:

$$m_{d_{BKE}}=2.75log({\tau }^{1.2}·0.37)-2.35.$$

(2)

This station, located in proximity of the crater axis, i.e., closer to where most earthquakes occur, was significantly more sensitive to the lower magnitude events. Indeed, where they overlap, the BKE catalogue included many more earthquakes with respect to OVO (Figure 3). As a consequence, this station was usually considered at INGV-OV to estimate $m_d$ for earthquakes recorded at the local seismic network and this was the value reported in the ALL catalogue.

Apart from the difference between these two catalogues, the number of events per year in the OVO catalogue significantly dropped after 2000, from an average of 318 earthquakes per year in 1972–2000 to 114 earthquakes per year in 2001–2021. In order to verify if this was a natural effect or if it was due to an oversight of the OVO catalogue, after the installation of the more sensitive BKE station in 1999, we compared the databases by excluding the earthquakes with $m<2.0$ (Figure 4), a value reasonably above the magnitude of completeness for both, as determined by the following analysis (Section 3).

The number of earthquakes per year in the two catalogues was very similar, differing by one at most, indicating that all of the local earthquakes with $m\ge 2.0$ for 1999–2021 were included in OVO; hence, this was likely to also be true for the preceding years. It follows that, rather than an artefact, the apparent decrease in the number of reported earthquakes after 2001 must have been a real effect, associated with reduced seismic activity within the volcanic structure. Thus, we used the whole catalogue for the statistical analyses. Besides these relatively long-standing instrumental earthquake catalogues, we also considered the ordinary, current catalogue (ALL; up to 9 October 2023), composed of earthquakes recorded from the whole monitoring network. It includes both located earthquakes and events for which no location could be determined. Thus, pointing to a full exploitation of the information, we also selected a subset of this database by extracting only the located earthquakes (LOC catalogue) and we analysed them separately. Because of the evolution of the network in the last 25 years (Figure 2), the number of located earthquakes increased significantly with time, in particular starting from 2015. This implied an improvement in the network sensitivity, with a decrease in the completeness magnitude.

3. The Gutenberg–Richter Distribution and the b-Value: Evidence for a Double Scaling

As a first step, we evaluated the completeness magnitude $m_c$ of the whole catalogues.

Several methods have been proposed for estimating this parameter. For instance, commonly used techniques are (i) the maximum curvature technique [75], which fixes $m_c$ at the magnitude value for which the maximum of the non cumulative GR distribution is observed; (ii) the goodness of fit [75], based on the exponential fit of the GR as a function of a threshold magnitude $m_{th}$, which selects $m_c$ as the first $m_{th}$ for which an accurate fit is obtained; (iii) the b-value stability approach [76], which chooses $m_c$ as the magnitude threshold value at which the estimated b value becomes stable; (iv) the entire magnitude range method [77], which uses the entire range of magnitude and multiplies the GR distribution by a complementary error function; and (v) the harmonic mean method introduced by [78,79].

Here, we adopt the method based on the evaluation of the variability coefficient $c_V$, which is defined as the ratio between the standard deviation of the magnitude and its average value [80] with $c_V$ = 1 for an exponential distribution. By evaluating its value as a function of a threshold magnitude $m_{th}$, it is possible to identify $m_c$ as the $m_{th}$ value, where $c_V$ reached a provided value. Here, we fixed this threshold value at 0.93, which simulations with synthetic catalogues revealed to be the best choice for retrieving the correct $m_c$, while smaller and larger values tended to underestimate and overestimate $m_c$, respectively. Figure 5 shows $c_V$ as a function of $m_{th}$ for the four analysed catalogues. The results indicate $m_c=-0.2$ for the BKE and ALL catalogues and $\sim 1.4$ for OVO (

Table 2. The completeness magnitudes, the b-values with their standard deviations, and the number of events N after the removal of the earthquakes with $m<m_c$.

Catalogue	${\mathit{m}}_{\mathit{c}}$	b	${\mathit{\sigma }}_{\mathit{b}}$	N
BKE	2.3	1.80	$0.15$	131
OVO	2.3	1.87	$0.76\times {10}^{-1}$	548
ALL	2.3	1.86	$0.75\times {10}^{-1}$	548
LOC	2.3	1.79	$0.14$	136
BKE	−0.2	0.75	$0.52\times {10}^{-2}$	15,631
OVO	0.9	0.79	$0.61\times {10}^{-2}$	7565
LOC	0.0	0.68	$0.68\times {10}^{-2}$	6895
ALL	−0.2	0.76	$0.55\times {10}^{-2}$	24,498

). Conversely, for LOC, $c_V$ overcomes the threshold value at $m_{th}\sim 1.8$, even though at $\sim 0.0$, it reaches a stable value. However, the break in the curve, with the exception of the OVO catalogue, reveals a more complex situation. Indeed, for three catalogues out of four, the GR distribution exhibited a two scaling pattern (Figure 6) confirming the results of [81]. As for OVO, we speculate that the higher $m_c$ might have possibly hidden the double scaling, which was clearly visible instead when the seismic moment $M_0$ (thus a non-logarithmic scale) was considered (see below, Section 3.1).

It is worth noting that, as expected, the GR distributions for ALL and BKE appearred to be almost identical (Figure 6). This was derived from the high BKE sensitivity to the lower magnitude Mt. Vesuvius seismicity, resulting in practically all of the occurring local earthquakes being included in the BKE catalogue. Indeed, the difference was mostly due to the ∼1400 events occurring during 2022–2023, when they did not overlap.

To better display the existence of the two regimes in the GR distributions and to estimate the two $m_c$ values, we evaluated $c_V$ in narrower magnitude ranges with respect to what is shown in Figure 5. The first range went from −0.5 to 1.8 and the $c_V$ evaluation confirmed $m_c=-0.2$ for the BKE and ALL catalogues (Figure 7). Conversely, OVO and LOC never reached the $c_V$ threshold value of 0.93. However, the $c_V$ value assumed an approximately constant value at $m_{th}=0.9$ for OVO and at $m_{th}=0$ for LOC. These values can be assumed to correspond to $m_c$ for those catalogues. The analysis in the second range, from 1.7 to 2.5, revealed that $m_c=2.2$ was true for all of the catalogues (Figure 8).

The analysis of $c_V$, in the two different magnitude intervals confirmed the existence of the two scaling regimes with different b-values. This feature appears to be very stable in time. Indeed, the GR distributions for each year (Figure 9) exhibited a double-scaling behaviour (although not very pronounced in some cases), which called for a deeper investigation on the possible causes of the experimental observation.

3.1. The Double Scaling of the b-Value: A Physical Effect?

The b-value was strongly affected by the properties of the crustal structure; in particular, stress conditions and fluid circulation, which also contributed to determining the patterns of the stress seismic release, in time and space. Thus, in the search for the origin, we first assumed that the observed double scaling in the b-value could result from the physical characteristics of the seismogenic volume. Potential sources could be, for instance, a finite size effect or the concurrence of earthquakes occurring in different sectors of the crustal structure, characterized by distinct properties, as also suggested by some authors [81].

As a first possible interpretation of the two scaling regimes characterizing the GR distribution, we checked if this could be attributed to a finite-size effect. Namely, the second (higher b) regime could be due to a tapering of the distribution with a faster decay, due to the system reaching its size limit. Following this hypothesis, the seismic energy should follow a Gamma distribution [82,83]. Thus, we transformed the magnitude to the seismic moment (proportional to energy), through the relationship $log{M}_0=9.8+0.8m_d$ [84]. The existence of two scaling regimes was also confirmed for the $M_0$ distribution for all of the catalogues (Figure 10), including OVO, revealing that, although hidden by incompleteness problems, the first scaling regime characterized this catalogue too. The existence of two different power law regimes in $p\left(M_0\right)$ led us to exclude the finite size effect as a possible explanation for the two scaling regimes in the GR distribution.

Then, we considered the hypothesis that the two distinct b-values could be associated with different earthquake clusters, separated in depth. In fact, two different physical mechanisms were proposed to cause shallow and deep events at Mt. Vesuvius, by indicating the gravitational load within the volcano edifice and a concurrence of regional stress and crustal heterogeneity (in association with the action of hydrothermal fluids) as the main sources of stress for shallow (depth < 1 km) and deep earthquakes [81], respectively.

We verified this possibility by evaluating the GR distribution for subsets of events in the LOC catalogue, either by repeatedly splitting the data according to a depth threshold $z_{th}$, varying by 0.5 km in a range of 0.5–2.0 km, or by separating the events into 1 km-thick layers. The results indicate that the deeper seismicity assumed smaller b-values in all cases (Figure 11), in agreement with previous studies [81] and with the evidence of a larger stress drop for deeper events [84]. The maximum difference in the b-value results by splitting the data at the depth threshold of 1 km. In fact, this level marked the separation between two major clusters of earthquakes —with the shallow one corresponding to events occurring inside the volcanic cone—possibly connected with distinct triggering conditions (e.g., [81,84]).

However, we remark that the GR distribution was characterized by double scaling for all of the $z_{th}$ values (Figure 12 and

Table 3. The number of events in the shallower and deeper layers n, the b-values with their standard deviations for the different $z_{th}$ values.

	Up			Down
${\mathit{z}}_{\mathit{th}}$	$\mathit{n}$	$\mathit{b}$	${\mathit{\sigma }}_{\mathit{b}}$	$\mathit{n}$	$\mathit{b}$	${\mathit{\sigma }}_{\mathit{b}}$
0.5	6086	0.79	$9.51\times {10}^{-3}$	2711	0.45	$5.94\times {10}^{-3}$
1.0	6901	0.75	$8.45\times {10}^{-3}$	1896	0.41	$5.91\times {10}^{-3}$
1.5	7434	0.71	$7.63\times {10}^{-3}$	1363	0.40	$6.67\times {10}^{-3}$
2.0	8038	0.67	$6.74\times {10}^{-3}$	759	0.41	$9.61\times {10}^{-3}$

) and in all the layers (Figure 13), leading us to exclude a depth dependence of b as a viable explanation for the observed break in the cumulative earthquake frequency-magnitude distribution.

Then, we investigated the dependence of the magnitude distribution on the distance from the Vesuvius crater (notably, the crater coordinates were extremely close to the average latitude and longitude of the located earthquakes). We separated the entire seismicity into two groups: one including the events inside a circle with a radius of r, centred at the crater, and the other with the remaining ones. As the threshold radius, we tested either the average ($d_a$) or the median ($d_m$) of the earthquake distance from the crater. In both cases, the GR law exhibited the same double scaling (Figure 14).

As a further tentative interpretation, we investigate the possibility that the two-scaling regime could result from the superposition of earthquakes characterized by distinct time clustering. In general, volcano-tectonic earthquakes occur in swarms—i.e., group of small earthquakes concentrated in time and space—or as variable mainshock–aftershocks sequences. Several authors considered that the different characteristics of these two groups might have resulted in distinct b-values (e.g., [85]). Therefore, we tested if this was the case for the seismicity of Mt. Vesuvius.

We verified this occurrence by separating the ALL catalogue in events occurring within and outside seismic swarms, respectively. Here, we focused only on the ALL catalogue, which was very similar to BKE, while OVO was almost Poissonian (with only two swarms individuated) and LOC was too scanty for this analysis. The method for individuating the swarms is described in the Appendix A.

The results illustrated in Figure 15 reveal that the GR distribution for both groups of earthquakes was very similar and, again, they were still characterized by the same two-scaling regime.

Apparently, the break of the slope observed in the frequency-magnitude of the Mt. Vesuvius seismicity was persistent with time and could not be explained as the superposition of earthquakes occurring either in distinct volumes or with different time clustering properties. Thus, based on the above analyses, we rejected the possibility that the double scaling in the b-value could be a real effect.

3.2. The Double Scaling of the b-Value: A Magnitude Uncertainty Effect

The use of the GR distribution for seismic hazard assessment or for deriving indications on the state of stress counted on the reliable determination of the earthquakes’ magnitude. For instance, Uchide and Imanishi [86] demonstrated that underestimation of the magnitude for microseismicity influenced earthquake statistics, thus affecting the seismic hazard inferences.

As described above, a duration magnitude $m_d$ for the Mt. Vesuvius seismicity was routinely estimated at the two sites of OVO and BKE, based on a relation derived in early 1980s from earthquakes with m ≥ 3.0, then extrapolated to lower magnitudes, typical of most earthquakes in this area.

Since then, the local seismic monitoring network run by the Osservatorio Vesuviano, now INGV, has been developed in the number and quality of the instruments for detecting, measuring, and recording the ground motion. Since 1998, Helicorder drums were progressively substituted with video monitors, while the number of station progressively increased since 2004, from 6 to the present 18 (Figure 2). The limited magnitude range for which a reliable duration–magnitude regression could be obtained and the technological changes might have had significant effects on the estimated $m_d$, possibly mirroring the resulting b-value.

For instance, Iannaccone and co-authors [87] demonstrated that the duration was overestimated by a few seconds, increasing with the waveform duration, when the seismograms were analysed on screen, rather than on drums. Because of the logarithmic dependence on the duration, this was reflected in a larger $m_d$—the difference exponentially increased with the decreased magnitude—with $\mathsf{\Delta }$$m_d$ being in excess of 0.1 and 0.2 at $m_d$ = 2.0 and $m_d$ = 1.0, respectively. In principle, being larger for small events, such a difference could determine the presence of a misleading double scaling in the b-value. However, if this was the source, it would require that the duration of a large majority of the events included in the analysed catalogues were estimated on Helicorder drums, whereas the complete transition to video monitor analysis of the seismograms was accomplished by 2000 [88].

Therefore, a different origin needed to be searched. To this aim, we focussed on the analysis of the regressions themselves. As mentioned above, the first duration–magnitude relation for Mt. Vesuvius seismicity was obtained for the OVO station, by comparing the observed duration of the ground motion associated with earthquakes occurring during the 1980–1981 Irpinia seismic sequence—50 to 80 km to E-SE—with the Richter magnitude $m_L$ estimated for the same even from the waveform recorded at a Wood–Anderson seismometer installed about 200 km NE of Vesuvius. The relative location of the sources and the recording sites did not allow for the calibration of earthquakes with a magnitude lower than 3.0; the regression was simply extrapolated to a lower magnitude [74], a procedure that could not ensure the correctness of the estimation at this range. As a matter of fact, Del Pezzo and Petrosino [89] demonstrated that this relation was not appropriate for smaller earthquakes. By using permanent and temporary stations, those authors derived the following relation:

$$m_L=0.682+0.655·m_d,$$

(3)

between the OVO duration magnitude and the BKE local magnitude, for 131 Vesuvius earthquakes, this showed that $m_d$ significantly underestimated $m_L$ at lower magnitudes (Figure 16). According to this relation, the duration magnitude at OVO should be evaluated using [89]:

$$m_d=1.8·log\tau -0.9.$$

(4)

We remark that the results obtained by [89] were well constrained, and were based (i) on the determination of an experimental distance–attenuation curve for Mt. Vesuvius events and (ii) on the simulation of Wood–Anderson seismograms, starting from recordings at a 1 Hz MARK L4-3D installed at BKE. On the other hand, less than 10% of the data analysed by [89] were above $m_d>2.0$, and it is quite apparent that, beyond this threshold, a simple $m_L=m_d$ regression (i.e., the one obtained by [74]) clearly justified the observations much better than the above relation.

Based on this observation, we split the whole dataset at $m=2$ and separately fit a linear model to the two groups of data. The regression for $m\ge 2$ resulted very close to $m_L=m_d$, being:

$$m_L=-0.18+1.09·m_d.$$

(5)

Conversely, the observation for $m\le 2$ are described by the following:

$$m_L=0.74+0.57·m_d,$$

(6)

very close to Equation (3), as shown in Figure 16.

The similarity of Equation (1) to Equation (5) implies that the relation derived by [74] for $m_d\ge 3.0$ could be reliably extended down $\sim 2.0$. At the same time, the difference using Equation (6) evidenced its fallacy for lower magnitudes. Thus, more reliable magnitude estimates for Mt. Vesuvius earthquakes could be obtained by using the relation (1) for $m_d\ge 2.0$ and the regression (6) for a lower magnitude.

On these grounds, we now focussed on the ALL catalogue and modified the magnitude by accounting for the above correction; then, we re-computed the frequency-magnitude distribution.

The corrected catalogue exhibited a more usual GR distribution (Figure 17), with a single, well defined, $m_c$ and a higher b-value, as generally observed in volcanic areas [90,91,92]. All of the above evidence strongly suggest that the break in magnitude scaling was not a real effect, rather it resulted from an arbitrary, inappropriate, extension of the regressions (1) and (2) to lower magnitudes. In order to retain the homogeneity of the catalogue, a change in the relation routinely used for the magnitude estimation is not recommended; however, we emphasize that caution should be provided when the seismic catalogues available for Mt. Vesuvius are used for investigating the time and space variations in the b-value. The relatively scarce level of the local seismicity would force the use of lower magnitude events in this kind of analysis, but their use appeared to not be appropriate below $m\sim 2.0$, if the magnitude was assumed as it was in the catalogues.

4. Time Variation of the b-Value and Temporal Clustering Properties

Once we correct for the bias in the magnitude estimation of the low energy events, we reevaluated the GR distribution throughout the considered time period and determined the time variations in the b value: an important tool for investigating the stress state changes with time. In addition, we also analysed the time variations of $m_c$, the average time interval between two successive events $\langle \mathsf{\Delta }t\rangle$, and the $\mathsf{\Delta }t$ variability coefficient $c_V^{\left(\mathsf{\Delta }t\right)}$, defined as

$$c_V^{\left(\mathsf{\Delta }t\right)}=\frac{{\sigma }_{\mathsf{\Delta }t}}{\langle \mathsf{\Delta }t\rangle },$$

(7)

where ${\sigma }_{\mathsf{\Delta }t}$ is the $\mathsf{\Delta }t$ standard deviation. This $c_V^{\left(\mathsf{\Delta }t\right)}$ should not be confused with the $c_V$ previously used in estimating $m_c$. Indeed, that is the magnitude variability coefficient, while here the $c_V^{\left(\mathsf{\Delta }t\right)}$ of $\mathsf{\Delta }t$ describes the time clustering properties of the seismicity; it assumes a value of 1 when the earthquake occurrence is Poissonian and larger values when events cluster in time. It is worth recalling that the $c_V^{\left(\mathsf{\Delta }t\right)}$ depends on the time of occurrence of earthquakes; thus, in principle, its analysis might be influenced by the assumed magnitude definition, affecting the $m_c$ value.

Figure 18 shows the four parameters as a function of time. $m_c$ was almost fixed at the two values of 0.6 and 0.7, in very good agreement with what was found for the whole catalogue. The b value exhibited a slight increase with time, indicating a possible decrease in stress on the volcano edifice. After removing the trend, b fluctuated around 1.4, confirming the results obtained for the whole catalogue. The b distribution was very close to a Gaussian one (Figure 19), whose standard deviation (0.13) was close to the average estimation error of b (0.08). This ensured that these fluctuations were purely random and not linked to the volcano dynamics. Also, $\langle \mathsf{\Delta }t\rangle$ exhibited few significant variations in fluctuation, within the standard deviation, around ≃18 h. On the other hand, $c_V^{\left(\mathsf{\Delta }t\right)}$ varied between 0.9 and 1.9; however, a great part of the $c_V^{\left(\mathsf{\Delta }t\right)}$ values fluctuated in the range of [1.1, 1.4], indicating a weak clustering degree (see the average rate of occurrence in the Appendix A). $c_V^{\left(\mathsf{\Delta }t\right)}$ assumed larger values corresponding with larger swarms.

The stationary behaviour of the analysed statistical parameters is a clear indication that the volcano state had not changed during the last 40 years. The only exception was represented by the occurrence of small swarms rapidly expiring (see Appendix A).

5. Conclusions

We analysed the seismicity of the Vesuvius volcano investigating the statistical characteristics of three different catalogues. The GR distribution appeared to be characterized by the presence of two scaling regimes unusual for the standard GR distribution. Here, we tested four different possible origins for this effect: (i) a finite size effect, (ii) depth variations of the b-value, (iii) a b-value depending on the distance from the crater area, and (iv) different b-value for time clustered and non-clustered events. None of these hypotheses are able to account for the double scaling regime of the GR distribution. As a consequence, we maintained that there was a bias in the magnitude. In particular, we suggest that the adopted duration–magnitude linear regression at the base of the catalogue, and originally determined for $m\ge 3.0$, cannot be extrapolated to smaller magnitudes. Thus, we derived a different relation between $m_L$ and $m_d$, for earthquakes with $m\le 2.0$, by using the data estimates by [89], and we corrected the catalogue magnitudes accordingly. When the correction was applied, the GR distribution assumed the usual shape and the resulting b was similar to the values usually observed in volcanic areas [90,91,92].

We also investigated the time variation of four statistical parameters: b-value, $m_c$, $\langle \mathsf{\Delta }t\rangle$, and $c_V^{\left(\mathsf{\Delta }t\right)}$. Our analysis indicates that the magnitude of completeness $m_c$ was not strongly affected by the evolution of the seismic network. This result was associated with the very localized nature of the seismicity around the crater axis, allowing for a good level of detection, even with the first, sparse network configurations. Similarly, the average time interval between successive earthquakes $\langle \mathsf{\Delta }t\rangle$ exhibited an approximately stationary behaviour, revealing that the time occurrence of seismicity at Mt. Vesuvius had not changed significantly in the last decades. The small changes of $c_V^{\left(\mathsf{\Delta }t\right)}$, excluding the abrupt increase in the correspondence of the three largest swarms, confirmed this result.

The analysis of the occurrence rate of the events within the swarms evidenced a sharp decrease with time (Figure A2), indicating the overall short duration for most of the almost 400 swarms, generally characterized by 10–20 earthquakes, with only 21 having more than 50 events.

Finally, the average Gutenberg–Richter b-values slightly increased during the entire analysed period, suggesting a slow, but continuous reduction in differential stress within the volcano edifice, possibly following the 1999 swarm, associated with the strongest earthquake (9 October 1999; $m_d$ = 3.6) since the last eruption in 1944–1945.