Southern-Tyrrhenian seismicity in space-time-magnitude domain

An analysis is conducted on a catalogue containing more than 2000 seismic events occurred in the southern Tyrrhenian Sea between January 1988 and October 2002, of minimum observed local magnitude 1.5, as an attempt to characterise the main seismogenetic processes active in the area in space, time and magnitude domain by means of the parameters of phenomenological laws. We chose to adopt simple phenomenological models, since the narrow space-time area did not allow to use more complex laws. The two main seismogenetic volumes present in the area were considered for the purpose of this work. The ﬁrst includes a nearly homogeneous distribution of hypocentres in a NW steeply dipping layer as far as about 400 km depth. This is probably the seismological expression of the Ionian lithospheric slab subducting beneath the Calabrian Arc. The second contains hypocentres concentrated about a sub-horizontal plane lying at an average depth of about 10 km. It is characterised by a background seismicity spread all over the area and by clusters of events that generally show a direction of maximum elongation. The estimated parameters of the statistical models describing seismogenetically homogeneous subsets of the earthquake catalogue in the three analysis domains, along with their conﬁdence intervals, are reported and analysed to establish whether they can be regarded as representative of a particular subset.


Introduction
The main task of this paper is the attempt to discriminate the seismicity relative to the different seismogenetic processes active in the Southern-Tyrrhenian area and characterise such processes by the parameters of mathematical models defined in the space, time and magnitude domains.This study requires the most important processes regulating the seismic release in the area to be preliminarily recognised.
The area (fig. 1) in which the epicentres are located includes the hinge zone between the emerged portion of the Sicilian Maghrebian Chain and the part of Southern-Tyrrhenian Basin featuring the highest crustal thinning.It lies within a wide W-E oriented shear zone cut by neotectonic NW-SE and NE-SW fault systems (Giunta et al., 2000).
The seismicity of the Southern-Tyrrhenian Basin is marked by a high frequency of events, average local magnitude about 2.2 and energy release rate about 2⋅10 12 J/yr.
The study was conducted on 2131 earthquakes (fig.2), each recorded by, on average, 8.9 stations of the Italian national seismic network and of some local networks, from January 1988 to October 2002.The hypocentres were determined with a procedure based on the joint optimisation of mean station residuals, hypo-central parameters and velocity model (Capizzi et al., 2001).The 481 epicentres marked in red in fig. 2 represent relatively deep events, which are clustered about a N-W dipping regression plane with an inclination of about 60° (Giunta et al., 2004).They are probably connected with the Ionian lithospheric slab which subducts beneath the Calabrian Arc.The black dots indicate the epicentres of the 1650 shallow events clustered about a sub-horizontal regression plane, on aver-age about 10 km deep, plunging both northwards and eastwards inside the thinned Southern-Tyrrhenian crust.Therefore, they are strictly connected with the strain field acting in the already mentioned hinge zone.
The detected seismicity of the Southern-Tyrrhenian hinge zone consists mainly of aftershock sequences, rarely preceded by foreshocks.Rather intense is also a background seismicity brought about by isolated events.The tendency of the shallow Tyrrhenian earthquakes to form clusters was gauged through an analysis of the correlation dimension of epicentres inter-distance (d2s) and inter-time distributions (d2t).The experimental information about hypocentral depth was not included in this analysis since the larger uncertainty in the estimation of this parameter would have produced an increase in the correlation dimension variance.The values of d2s, determined with the Maximum-Likelihood (ML) and the Least-Squares estimators (LS) defined and analysed in De Luca et al. (1999Luca et al. ( , 2002)), are 0.50 and  (Giunta et al., 2004).( Giunta et al., 2004).

Clusters and independent events
To detect each single cluster of strongly correlated events and the isolated ones in a seismic catalogue, an algorithm was devised based on the assumption that a dependent event can be generated if the stress-field perturbation produced by a previous event goes over a given intensity threshold.Making the simplistic assumption that the perturbation intensity monotonically decreases with its space-time distance from the generating event, it appears reasonable to hypothesise that the probability of a subsequent event to be dependent on the previous one decreases as well.
Choosing a threshold beyond which two events are to be considered belonging to different clusters is equivalent to defining two thresholds values δ S and δt for the space and time distances from the former event, respectively.These parameters ought to be defined according to the event magnitude and to the seismogenetic volume rheology.Even if such choice could depend on the magnitudes, because of the limited observed range in the space, time and magnitude domains, two single threshold values for all magnitudes are chosen and a simple windows-based clustering method (see for example Reasenberg, 1985) is used.
In the analysis of a catalogue, therefore, an event i is attributed to a cluster if at least another event j is found such that the double condition and holds, where ∆ S(i, j) and ∆ t(i, j) are the interdistance and inter-time between the two events, respectively.
The two threshold values were chosen looking at the distribution of the inter-points distances, also taking into account the main geological structures in the area, and were fixed to δS=35 km and δt =2 days.The extraction procedure leads to the determination of 185 clusters, among which 8, represented in the map in fig.3, contain more than 10 events each; their space distribution suggests a possible correlation with the main fault systems reconstructed on the basis of the geological information and the con- # δ ∆ 0.55 respectively, whilst d2t, computed by the least-squares method, is 0.74.These estimates, being both significantly lower than the embedding dimensions in their respective domains (2 and 1), indicate an evident tendency of the Tyrrhenian seismogenesis to gradually release strain energy.
The distribution of faults, stress field, geometry of seismogenetic volumes and their depth, this latter affecting pressure and temperature, control the way in which strain energy is released, through clusters and isolated events.For this reason one might expect to observe significant differences in the parameter estimates of a phenomenological model applied to sets of earthquakes relative to different release modes.Obviously, the significance of these differences should be assessed by statistical tests.
Such an analysis requires that the seismic catalogue is first divided into subsets of clustered or isolated, shallow or deep events.The clusters and the background events, both present in the same seismogenetic volume, can be separated either by defining a priori a complete set of differentiation criteria or choosing a multi-domain objective function (e.g., a likelihood function), which can be optimised by iteratively changing the partitioning of the events (Adelfio et al., 2006).
Therefore, in the next section the windowsbased clustering technique used is described along with some descriptive results relative to the obtained partition; in the third section the analysis in the space domain is carried out, mainly focussing on a non-parametric method to estimate the spatial density of events.The fourth section is about the analysis in the magnitude domain, with some empirical analysis for the determination of the magnitude completeness threshold and an attempt at comparison between some of the homogenous sets of earthquakes disclosed by the clustering method used; the time domain analysis is performed in the fifth section both for the main aftershocks sequences (through the estimation of the Omori's laws parameters) and for the set of isolated events (with the study of stationarity in time); finally in the sixth section we report final remarks and some issues for future researches.straints set by the kinematic model assumed for the area.
In particular, the blue cluster in fig. 3 contains 541 earthquakes and is related to the seismic sequence that followed Palermo's event which occurred on 6 September, 2002; the green one contains 69 earthquakes which occurred in the offshore of San Vito Lo Capo (Trapani) in June 1998.These two clusters alone have a numerosity that allows a significant statistical analysis.
Moreover, an intense background seismicity, constituted by the isolated events and by the main shock of each sequence, is present in this zone.
A description of the partition obtained is provided by the estimate of the space and time correlation dimension of the set of independent events and of the complete set.The values of d2s, determined by the ML and LS estimators for the independent events (De Luca et al., 1999, 2002), are on average 1.59 and that of d2t is 0.97, whilst the corresponding ones for the complementary set (i.e. the set of the clustered events) are 0.42 and 0.75, respectively.To discriminate the two seismicity components in the space-time domain, defined here as the Cartesian product of the two spatial dimensions and the temporal one, a parameter indicating the total clustering degree d2st can be defined as the sum d2 s+ d2t, which turns out to be 2.56 and 1.17 for the two sets.This parameter suggests a sharp difference between the phenomena originating the two extracted components.
The fractal characterisation in space and time domain of single aftershock clusters does not differ significantly from that of independent events.As an example, the d2s obtained for the cluster with 541 events and the one with 69 events are on average 1.61 and 1.87, respectively, while the d2t are on average 0.89 and 0.73, corresponding to d2st equal to 2.50 and 2.60, respectively, not significantly different from the d2st of the independent component.This result is probably justified by the fact that the location errors and the low numerosity of observed clusters do not allow us to observe structures at a scale lower than that of the whole cluster.(Giunta et al., 2004).

Analysis in space domain
In order to define a number of geometric parameters useful to distinguish seismically homogeneous areas through the rate of seismic release relative to a given seismicity component (independent or cluster), it is useful to determine a continuous function expressing the intensity of the non-homogeneous point process.In particular, a non-parametric estimation technique is used in this work in order to estimate a density function on the basis of a sample P 1, P2, ..., PN of N epicentres (Xi,Yi), being unable to hypothesise a functional form f(x,y;θ) for the observations density.
The simplest and most widespread method for a density estimation consists in the construction of the frequency histogram of observations within equally sized rectangular bins.Even though this method is very useful for a preliminary description of the data, the estimate of a continuous density function is often preferable.This latter can be obtained through interpolation and smoothing procedures, amongst which, frequently adopted are those based on kernel estimators.The kernel estimator is defined as where K(x, y; X i, Yi, hx, hy) is the generic kernel function centred at the point (X i, Yi) and (hx, hy), called smoothing parameters, are constants expressing the degree of peakedness of kernel functions.The function f has to be constrained by the positivity and normalisation conditions in its domain, which allow to regard it as a density function; also, the continuity and differentiability properties of the function K must hold.The kernel estimator can be obtained as the sum of N surfaces, each corresponding to the single pair of epicentral coordinates: the kernel function defines the shape of these surfaces.
The application of this technique requires the parameters hx and hy to be optimised in order to achieve the best compromise according to some subjective criteria between the surface smoothness and the level of detail in the phenomenon representation.

/
To tackle the optimisation problem, one has to be aware that, as h x and hy grow, the estimator variance decreases but its bias increases and the density estimate within the investigated area reaches a slightly decreasing value (the estimator is progressively more affected by further observations).The optimum value of (hx, hy) is chosen such as to minimise the random mean square error of the estimator f ˆ(⋅) around the true density f(⋅), i.e.
The choice of the kernel should lead to a mean square error as low as possible, under the assumption that the smoothing parameters (h x, hy) have been chosen correctly.Since it is known that the efficiency of the density estimator (in terms of its MSE) poorly depends on the adopted kernel function K(⋅) (Silverman, 1986, p. 43), the kernel selection is made on the basis of different criteria, such as the requested degree of differentiability or the computational burden.The estimator that we adopted on the basis of such criteria has the bivariate Gaussian kernel as K-function.The smoothing constant is evaluated with Silverman's (1986) formula, which optimises the estimator asymptotic behaviour in terms of mean square error and provides valid results on a wide range of distributions.It is obtained from with A=min{standard deviation, range-interquartile/1.34}.
It turns out from the analysis that the values in km obtained for hx and hy are 0.80 and 1.10 for Palermo cluster, 2.7 and 2.6 for San Vito Lo Capo cluster and 21.7 and 7.1 for the independent events.The strong difference between hx and hy observed for the independent events can be attributed to a lower variability of the epicentral distribution along E-W than N-S.The small differences shown by the estimates of hx and hy for Palermo and San Vito Lo Capo clusters reflect the nearly NE-SW and NW-SE orientations of the fault systems that generated them., ). var . Estimated spatial density of the independent events.The density distribution for the set of independent data (fig.4) evinces a maximum intensity of seismicity in the zone surrounding Patti (which could however be biased by the higher network sensitivity in the area) as well as in a band sub-parallel to the Ustica-Aeolian Islands line.
The intensity functions for the events of Palermo sequence are calculated both including all events and subdividing them into three subsets corresponding to three different depth intervals containing each approximately the same number of events: hypocentres lying within 7 km depth, in the interval [7,10] km and deeper than 10 km (fig.5).
Comparing the distributions relative to the three depth intervals, it is reasonable to conclude that all hypocentres are concentrated inside a narrow sub-vertical seismogenetic volume having a crustal thickness and a length about 40 km.More precisely, the hypocentres concentrate, with an mean square distance equal to 2 km, about a regression plane with strike 234°dipping 83°northwest.
The significance of the density distribution estimated for San Vito Lo Capo cluster of 1998 is limited by data insufficiency; despite that, it is possible to obtain a description of the intensity of the spatial phenomenon, characterised by a NW-SE stretching direction in agreement with Castellammare del Golfo system fault (fig.6).
To attempt a description of the main space attributes of sets of earthquakes through synthetic parameters, the magnitude of the mean horizontal gradient of the estimated density functions is evaluated.It turns out that the ratios between the values of this parameter estimated for the two main clusters and that relative to background seismicity are 1372 (Palermo) and 182 (San Vito Lo Capo), respectively.This result shows that, beyond the fact that typical mean horizontal gradients for clusters are over 100 times as large as those related to background seismicity, strong differences can be found between the gradients of single clusters, mostly connected with the geometry of release volumes in addition to the total seismic energy released.

Analysis in magnitude domain
The set of independent events, that of the Ionian slab and the two clusters are also analysed in the magnitude domain through the estimate of b of Gutenberg-Richter's Law, in order to evaluate the possibility to discern heterogeneous seismogenetic conditions.To achieve estimates of b not biased by the catalogue incompleteness, an attempt is made to determine the magnitude completeness threshold (M 0) for each of these sets.To estimate M 0 using uniform criteria two different methods are used.The first consists in the examination of the differences ∆b of consecutive estimates of b for increasing magnitude completeness threshold: M0 is the point up to which such differences are approximately constant.
The result of this test provides narrow intervals for the estimated thresholds.As an example, fig.7 shows ∆bML/∆M0 as a function of M0 for the set of independent events, where ∆bML is the difference between two consecutive maximum likelihood estimates of the parameter b.A nearly constant ∆b ML /∆M0 up to M0=2.5 can be observed in fig.7, followed by an interval in which it becomes strongly variable assuming a nearly zero average.
The second tool compares the experimental frequencies and the frequencies expected from the Gutenberg-Richter's Law for magnitudes greater than M 0, by using the Kolmogorov-Smirnov test.Denoting by x (i) the i-th order statistic, the expression of the Kolmogorov-Smirnov statistic is which determines the largest difference between the theoretical cumulative distribution F(⋅) and the empirical one Femp(⋅).From the analysis of the different P-values for magnitudes greater than 2.3, it is observed that they never reach values lower than 0.5 as the magnitude grows.A global evaluation of the tests carried out suggests putting 2.6 as threshold value.The same value is fixed for the set of Palermo cluster, whilst one included between 2.3 and 2.5 is regarded as suitable for the events of the Ionian slab.The lower threshold magnitude found for the slab events can be explained by a better average coverage of the seismic network than in the case of the other two sets.
The parameter b is estimated both through Utsu's (1965) formula, b ˆML, and by means of the estimators by Tinti and Mulargia (1987) b ˆTM and by Bender (1983), b ˆB, not biased by the binning of magnitude data into frequency classes of width 0.1, in order to evaluate the effect of different estimators on each set of events.The estimates are reported in table I.
The confidence intervals of the ML estimates are found through the asymptotic distribution of the Likelihood Ratio statistic (LR).In the statistical inference theory the LR is defined as a test statistic, asymptotically distributed as a χ 2 random variable, obtained by the ratio between the maximum of the likelihood function under a fixed null hypothesis and the value of the likelihood function corresponding to the ML estimator.
For the three sets of events defined above a comparison among the estimated values of b is carried out, even if they are computed with different magnitude completeness thresholds.To assess the null hypothesis that the b's of the Gutenberg-Richter's laws fitted to the three sets of earthquakes are equal, for the maximum likelihood estimates a Likelihood Ratio test is used.The value of the test statistic is 4.58 with a Pvalue 0.10.
The decrease of b, usually observed as hypocentral depth increases (Wyss, 1973), is not highlighted in the analysis of this dataset.

Analysis in time domain
Event sets generated by different seismogenetic processes could behave differently in time domain.An analysis is conducted on the events that make up background seismicity with the purpose of assessing the hypothesis after which such events constitute a Poisson process without memory in a temporal sense; this latter characterises phenomena in such a way that the probability of a future event is constant and not dependent on its past.
The homogeneous Poisson process represents the basis upon which the theory of point processes is grounded, and many processes can be obtained as one of its generalisations.It represents the simplest stochastic mechanism for the generation of a map of points, and is used as a standard reference for a case of complete space randomisation.
If the phenomenon is regarded as a time sequence of events, the point process {T n, n≥1} in R + is called Poisson process of intensity λ>0 if and only if its related counting process {N(t), t≥0}, which describes the number of events or successes in [0, t], satisfies the three conditions: 1) For each s, t≥0, N(t+s)−N(t) is a random variable with a Poisson distribution of mean λs con k∈N.
2) {N(t), t≥ 0} has independent increments, i.e. the numbers of events occurred in two disjoint time intervals are independent of each other.
3) N(0)=0.It should be observed that, in the one-dimensional case, the variables t1=T1, t2=T2−T1, ..., tn= =T n−Tn−1, ... constitute a continuous process, indexed by n ∈N, called waiting-time process, which is defined by independent and equally distributed variables, with an exponential distribution with intensity λ.Thus, as an event occurs, the process restarts as a replication of the initial one, featuring the so-called absence of memory effect.
In order to assess the assumption of independence for the set of the events previously recognised as isolated by the clustering procedure, an analysis of the inter-time distribution is conducted.To such purpose, the Kolmogorov-Smirnov statistic is computed to verify the fit of an exponential distribution with the empirical one of the observed inter-times.As a comparison, the same test is repeated for the complete set of earthquakes.The test statistic values estimated for the two sets of events are 0.04, with P-value 0.24, and 0.42 with P-value 2.20 ⋅10 −16 , respectively, therefore there is evidence with a 5% error in favor of the exponential distribution for the inter-event times relative to background seismicity, whereas such an assumption is to be rejected for the complete data set.Figure 8 re- " , ports the histogram of the observed relative frequencies and the theoretical distribution, resulting from the distribution analysis of the intertimes between independent events.For this set the estimated value of λ is 0.127.Moreover, as shown in fig.9, the mean of the inter-times between isolated events (dt's) does not change significantly in time.Indeed the estimated coefficient of the least-squares line fitted to points is −0.001, and carrying out a t-test, to assess if this value is significatively different from zero, the value of the test statistic is 1.909 with a Pvalue greater than 0.05.Therefore the temporal properties of the independent events seem to be consistent with those of a homogeneous Poisson process with a constant intensity λ.
As already known, the time evolution of a seismic sequence behaves like a relaxation    Southern-Tyrrhenian seismicity in space-time-magnitude domain process generally described by Omori's Law (Utsu, 1961), It is believed that the parameters of this law are affected by certain physical and geometric features of the seismogenetic volume; e.g., p has been observed to decrease with the heterogeneity and complexity of this volume.We apply this law to the two main clusters present in the catalogue, in order to assess the model adaptability to these sets as well as evaluate the significance of the differences between the estimates of the parameters involved.
The law parameters are estimated with the maximum-likelihood method (Console, 2001) and the results relative to Palermo cluster reported in table II for different threshold magnitudes M 0.
The estimated p takes values between 0.75 and 0.81 for the different M0.Comparing their 95%-confidence intervals, determined by means of the asymptotic theory, it turns out that there is no significant dependency of the estimates on the threshold magnitude.
The analysis of the test statistic X 2 and even more the P-value estimates, with a 95%-confidence (table II), indicate a bad fitting of the theoretical law to the observed frequencies.
In fact, the energy release for Palermo sequence occurred in a rather complex way.Twenty-one days after the beginning of the sequence, an event with local magnitude 4.6, one unit lower than the main shock, seems to have been followed by its own aftershocks sequence (fig.10); the latter, overlapping with the sequence generated by the main event, produced a sensible increase of seismic activity, as already observed in several other cases.
Describing the complete time evolution of such sequences would require definitely more complex mathematical models, which could be safely applied if a longer time interval was available.As consequence, the parameters of Omori's Law are also estimated on the sequence truncated immediately before the event with magnitude 4.6.The parameters and confidence intervals estimates (table III and fig.11) are significantly different with respect to those obtained for the non-truncated sample and the test statistic X 2 and the P-value indicate a neat improvement of the theoretical law fitting to experimental frequencies.Also in this case the estimated parameter p and confidence inter-vals do not show a significant dependence on the threshold magnitude.
In the analysis of San Vito Lo Capo sequence, understandably deprived of its foreshocks, two magnitude thresholds alone are taken into account because of its low number of events.In this analysis the p's (table IV) turn out much smaller than those of Palermo sequence, and the parame-ter estimates are scarcely significant, as evinced by the test of hypothesis rather than from the confidence intervals, extremely unstable because of the sample insufficiency.

Conclusions
This paper deals with a statistical analysis conducted on the seismicity of the Southern-Tyrrhenian Sea, aimed at establishing whether the experimental information available allows to discern, through parameters relative to wellknown phenomenological laws, heterogeneous seismogenetic processes which have been previously identified in a qualitative way.
The analysis in space domain indicated that, given the large amount of information currently available on the studied area, it is possible to set up mathematical models that describe the intensity of the seismogenetic process in a certain area with a sufficient level of detail.
As concerns the seismic release modes, a strong tendency was observed of the seismicity in the Southern-Tyrrhenian Sea to take place through sequences of aftershocks.Because of the insufficient degree of completeness and accuracy of low-energy experimental information it was not possible to analyse the degree of complexity of the seismogenetic process relative to each sequence.For this reason, the sets of aftershocks were characterised with parameters describing average properties of the estimated density functions.
The analysis in the time and magnitude domains pointed out the need to formulate more general, possibly multi-domain, mathematical models, capable of predicting the process complexity.Since such a description would require a wider space-time area than that observed, we resorted to the most commonly adopted phenomenological models.
The quality of parameter estimates seems to have been limited by the model inadequacy, yet also by the small amount of observations.In fact, the comparison between point and interval estimates of the investigated parameters suggested that the differences between the estimates of b and those between the parameters of Omori's Law for distinct sets of earthquakes are not statistically significant.It also turned out that all such estimates strongly depend on the adopted estimator, and their quality is affected by the completeness analysis on the data.

Fig. 5 .
Fig. 5.Estimated spatial density of Palermo cluster relative to all events and to each depth interval.

Fig. 7 .
Fig. 7. Variation rate of parameter bML as a function of M0 for the set of independent events.

Fig. 9 .
Fig. 9. Inter-times versus time for the independent events and fitted linear trend.

Fig. 8 .
Fig. 8. Histogram of the inter-times for the independent events compared with the exponential density function.

Fig. 11 .
Fig. 11.Observed aftershocks frequency compared with Omori's Law for the truncated Palermo sequence.

Table II .
Parameters of Omori's Law estimated for different threshold magnitudes M0, confidence intervals for p and c, X 2 and P-value of the X 2 -test (Palermo sequence).

Table III .
Parameters of Omori's Law estimated for different threshold magnitudes M0, confidence intervals for p and c, X 2 and P-value of the X 2 -test (truncated Palermo sequence).

Table IV .
Parameters of Omori's Law estimated for different threshold magnitudes M0, X 2 and P-value of the X 2 -test (San Vito Lo Capo sequence).