Nonextensive analysis of crustal seismicity in Taiwan

Using the Taiwan Central Weather Bureau earthquake catalogue, the crustal seismicity of Taiwan was analyzed by means of a nonextensive approach. The time span of the analyzed catalogue is from 1 January 1990 to 30 November 2007, and only earthquakes with magnitude M≥2.0 were considered. Our findings reveal that the nonextensive statistics furnishes a very good prediction of the cumulative magnitude distribution of crustal seismicity in Taiwan, even if the aftershocks are removed, indicating that the approach is robust for clustered as well as declustered seismicity.


Introduction
The well known empirical Gutenberg-Richter (1944) relationship between the cumulative number of earthquakes (for magnitude larger than a threshold M th ) with M th explains statistically the seismicity occurrence, but was not related with general physical principles.Sotolongo-Costa and Posadas ( 2004) developed a general model for earthquake dynamics, and the Gutenberg-Richter law could be considered a particular case.In this model the mechanism of earthquake triggering is based not only on the slippage of fault planes and relative displacement due to the break of the asperities (De Rubeis et al., 1996), but is also caused by the fragments filling the space between fault planes and originated from the breakage of the tectonic plates.Using the non-extensive Tsallis formalism (Tsallis, 1988), a more realistic magnitude distribution was deduced, providing an excellent fit to the seismicities of several seismic regions (Silva et al., 2006;Vilar et al., 2007), and from regional seismic catalogues down to specific seismic areas (Telesca, 2010a, b).

Nonextensive earthquake model
An earthquake fragment-asperity interaction model was developed by Sotolongo-Costa and Posadas (2004) and, then revised by Silva et al. (2006), in which the earthquake triggering mechanism is given by the interaction between the irregularities of the fault planes with the fragments between them, originated by the local breakage of the tectonic plates, from which the faults are generated.In this model the relative position of fragments filling the space between the irregular fault planes can contribute to the hindering of their relative motion.Then, stress increases until a displacement of one of the asperities, due to the displacement of the hindering fragment, or even its breakage in the point of contact with the fragment leads to a relative displacement of the fault planes of the order of the size ρ of the hindering fragment, producing a liberation of energy ε (Sotolongo- Costa and Posadas, 2004).As large fragments are more difficult to release than small ones, this energy is assumed to be proportional to the volume of the fragment, so that the energy distribution of earthquakes generated by this mechanism can reflect the volumetric distribution of the fragments between plates (Silva et al., 2006;Vilar et al., 2007).The use of nonextensive statistics is a suited tool to describe the volumetric distribution function of the fragments.Applying the maximum entropy principle for the Tsallis entropy (Tsallis, 1988), the Tsallis entropy for our problem is given by where p(σ ) is the probability of finding a fragment of surface σ and q is a real number (Vilar et al., 2007).the ad hoc condition about the q-expectation value, σ q = σ q = ∞ 0 σ p q (σ )dσ .For q→1, the last condition becomes the definition of mean value.
Applying the standard method of conditional extremization of the entropy functional S q , we obtain the following expression for the fragment size distribution: which corresponds to the area distribution for the fragments of the fault planes (Silva et al., 2006;Sotolongo-Costa and Posadasa, 2004).
Assuming the energy scale ε∼r 3 , the proportionality between the released energy ε and r 3 becomes σ − σ q = ε a 2 3 , where σ scales with r 2 and a (the proportionality constant between ε and r 3 ) has dimensions of volumetric energy density.This scale is in full agreement with the standard theory of seismic moment scaling with rupture length (Lay and Wallace, 1995).Thus, using the transformation σ − σ q = ε a 2 3 , the following form for the energy distribution function of the earthquakes follows: with and . The probability of energy p(ε)=n(ε)/N, where n(ε) is the number of earthquakes of energy ε and N the total number of earthquakes.Therefore, the cumulative number of earthquakes can be calculated as the integral of Eq. (3): Where N (ε > ε th ) is the number of earthquakes with energy larger than the threshold ε th .Due to the relationship m = 1 3 log(ε), where m is the magnitude, the following formula for the distribution of the number N of earthquakes whose magnitude m is larger than the threshold M th normalized to the total number of events is obtained: (5) This not trivial result incorporates the characteristics of nonextensivity into the distribution of earthquakes by magnitude.The q-value is a quantitative measure of the length scale of the spatial interactions: q∼1 indicates short-ranged spatial correlations; as q increases, the physical state becomes much more unstable; high values of q mean that the fault planes are not in equilibrium and more earthquakes can be expected (Sotolongo-Costa and Posadas, 2002).

Data analysis
The Taiwanese seismicity from 1 January 1990 to 30 November 2007 was investigated.The data were released from the Central Weather Bureau (CWB) earthquake catalogue (Fig. 1).The minimum magnitude M for which the catalogue can be considered complete is 2.0 (Telesca et al., 2009).
Figure 2 shows the normalized cumulative magnitude distribution (NCMD) of the Taiwanese crustal seismicity.The fitting of NCMD with Eq. ( 5) was performed by means of the Levenberg-Marquadt nonlinear least square method (Levenberg, 1944;Marquardt, 1963).When using this method, the initial values q 0 and a 0 are firstly fixed, then the algorithm furnishes the values for q and a that best fit the data and lie within admissible ranges that, in this study, were [1, 2] for q and [10 −3 , 10 12 ] for a.The best fitting is established after the minimum least square fitting error e has been reached.error e a 0 q 0 =1.1 q 0 =1.2 q 0 =1.3 q 0 =1.4 q 0 =1.5 q 0 =1.6 q 0 =1.7 q 0 =1.8 q 0 =1.9 254 Fig. 3(a Since a priori the initial values q 0 and a 0 giving the minimum estimation error e are not known, the error e was calculated varying the initial values of a 0 and q 0 .The error e is given by the square-2 norm of the residual, which is the difference between the original data and the values evaluated by the fitting nonlinear function.Figure 3 shows the error e, the parameters q and a, respectively, varying the values of q 0 and a 0 .The value of the error e is stable for a 0 ≥10 4 ; in this range the obtained values of q and a are also quite stable, with q around 1.68 and a of the order of 10 7 .The minimum error is ∼0.2868 and was obtained for q 0 =1.4 and a 0 =10 7 ; then, these initial values were used to find the best fitting parameters q=1.679 and a=1.4×10 7 .The fitting curve is also shown in Fig. 2 as a solid line.The agreement between the NCMD with Eq. ( 5) is very good, indicating that the nonextensive fragment-asperity interaction model is well suited to describe the crustal seismicity of Taiwan.
In the Taiwanese territory many strong earthquakes occur and the aftershock activity is, therefore, also quite intense.In order to verify if aftershocks could influence the values of q and a, the nonextensive analysis was also performed on the aftershock-depleted catalogue.For the algorithm of removing aftershock, we had applied the method of spatiotemporal double-link cluster analysis to eliminate aftershocks in the catalogue.This method is similar to the single-link cluster analysis proposed by Davis and Frohlich (1991).Given a magnitude threshold of main shocks, the declustering algorithm specifies two linking parameters in the time and space scales, 3 days and 5 km for example.An event would be identified as an aftershock when its epicenter and occurrence time lay within the prescribed spatiotemporal window of some main shock.Then, the procedure iteratively searches for secondary aftershocks, i.e. the aftershock of an earlier aftershock.The whole catalogue turns out to be separated by whole catalog error e a 0 q 0 =1.1 q 0 =1.2 q 0 =1.3 q 0 =1.4 q 0 =1.5 q 0 =1.6 q 0 =1.7 q 0 =1.8 q 0 =1.9 254 Fig. 3(a whole catalog q a 0 q 0 =1.1 q 0 =1.2 q 0 =1.3 q 0 =1.4 q 0 =1.5 q 0 =1.6 q 0 =1.7 q 0 =1.8 q 0 =1.9 q 0 =1.1 q 0 =1.2 q 0 =1.3 q 0 =1.4 q 0 =1.5 q 0 =1.6 q 0 =1.7 q 0 =1.8 q 0 =1.9 258 Fig. 3(c) 259 whole catalog q a 0 q 0 =1.1 q 0 =1.2 q 0 =1.3 q 0 =1.4 q 0 =1.5 q 0 =1.6 q 0 =1.7 q 0 =1.8 q 0 =1.9 q 0 =1.1 q 0 =1.2 q 0 =1.3 q 0 =1.4 q 0 =1.5 q 0 =1.6 q 0 =1.7 q 0 =1.8 q 0 =1.9 258 Fig. 3(c) 259 Fig. 3.The error e (a) and the parameters q (b) and a (c) versus initial value a 0 , calculated for the whole seismicity.The different colored curves are obtained with different values of q 0 .The value of the error e is stable for a 0 ≥10 4 ; in this range the obtained values of q and a are also quite stable, with q around 1.68 and a of the order of 10 7 .The minimum error is ∼0.2868 and was obtained for q 0 =1.4 and a 0 =10 7 ; then, these initial values were used to find the best fitting parameters q=1.679 and a=1.4×10 7 for the fitting curve, shown in Fig. 2. error e a 0 q 0 =1.1 q 0 =1.2 q 0 =1.3 q 0 =1.4 q 0 =1.5 q 0 =1.6 q 0 =1.7 q 0 =1.8 q 0 =1.9 many sequences of main shock and aftershocks.By using the temporal and spatial linking parameters of 3 days and 5 km, we had removed the aftershock events generated from main shocks with M L >4.5.Those linking parameters were usually used for declustering the CWB catalogue (Wu and Chen, 2007).
Figure 4 shows the NCMD for the declustered catalog along with its fitting curve, which was obtained similarly to the whole case.Figure 5 show the error e, q, and a varying the initial values of a 0 and q 0 , like in the whole catalogue.
The error e shows an almost flat behavior for a 0 ≥10 4 and for any q 0 value.This feature does not allow to estimate a minimum error e, in order to find the suited initial values of a 0 and q 0 , and, in turn, the value for the fitting parameters q and a.But, due to the very small variation of the obtained parameters q and a in the same range as the error e, it is reasonable to estimate these as the mean of the obtained values for a 0 ≥10 4 ; thus q=1.685 and a=1.1×10 7 ; such values are very close to those obtained for the whole catalogue, indicating that the nonextensive model is robust with the presence of aftershocks.

L
. Telesca and C.-C. Chen: Nonextensive analysis of crustal seismicity in Taiwan The maximum entropy formulation for Tsallis entropy implies that two conditions have to be introduced (Tsallis et al., 1998): 1) the normalization of p(σ ),

Fig. 1 .
Fig. 1.Spatial distribution of the Taiwanese seismicity from 1 January 1990 to 30 November 2007.The data were extracted from the Central Weather Bureau catalogue.

Fig. 5 .
Fig. 5.The error e (a) and the parameters q (b) and a (c) versus initial value a 0 , calculated for the aftershock-depleted seismicity.The different colored curves are obtained with different values of q 0 .The error e shows an almost flat behavior for a 0 ≥10 4 and for any q 0 value.The parameters q and a were estimated as the mean of the obtained values for a 0 ≥10 4 , thus q=1.685 and a=1.1×10 7 (see text for details).