Scaling Properties of Aftershock Sequences in Algeria-Morocco Region

the threshold completeness m c , we used the maximum likelihood to estimate above the threshold completeness m c . The results obtained are close to 1.0, the typical universal value for aftershock sequence. The Omori-Utsu law for aftershock decay and Bath's law for the difference between magnitude and the largest aftershock and mainshock as modified by Shcherbakov et al., (2004a, b)., while using AIC as a measure to select the most appropriate model between the first stage and two stage Omori-Utsu model, the rate of decay of aftershocks was found to follow in all case, except the May 21, 2003 Zemouri aftershock sequence, the first stage Omori-Utsu law, which mean the Omori-Utsu model without secondary aftershock. The decay of aftershock trigged by the May 21, 2003 Zemouri earthquake exhibit a better fit with a two stage Omori-Utsu model, denoted model 4 in the text. The later allow us to include in the modelisation the secondary aftershock of magnitude 5.8 m bLg . The study of the Omori-Utsu law has been complemented of aftershock sequences. This procedure is an alternative to the Omotri-Utsu law based on the physics of the static fatigue mechanism. The study of Guten‐ berg-Richter relationship and the modified Bath's law, allowed us to perform the partitioning of the energy released during each aftershock sequence by the main shock and aftershocks. The spatial analysis has been performed using the fractal dimension D 2 derived using the Correlation integrale. Using non parametric estimation approach, especially the kernel density estimation, we derive the estimation of the density of the probability of the inter-event distances. The mode of this density highlight how the spatial clusters of the events. This study is a first attemp to perform analysis of aftershock sequneces triggered by main events in the Algeria-Morocco region.


Introduction
This chapiter is dedicated to the analysis of some aftershock sequences occurred in Algeria-Morocco region, namely Al Hoceima earthquakes of May 26, 1994  Aftershock sequence is usually attributed to the strain energy not released by the mainshock. Statistical properties of aftershocks have been extensively studied for long time. Most of them dealt with the distribution of aftershock in time, space and magnitude domain. Several authors have noted the importnace of systematic investigation of aftershock sequences to earthquake prediction and a number of statistical models have been proposed to describe seismicity characters in time, space and magnitude (Utsu, 1961;Utsu et al., 1995;Guo and Ogata, 1997;Ogata, 1999). It is also well known that aftershock sequence offer a rich source of information about Earth's crust and can provide and understanding of the mechanism of earthquakes, because tens of thousands of earthquakes can occur over a short period in small area. The tectonic setting and the mode of faulting are factors others than fault surface properties that might control the behavior of the sequence (Kisslinger and Jones, 1991;Tahir et al., 2012). It is widely accepted that aftershock sequence presents unique opportunity to study the physics of earthquakes. In the same time importnat questions concern the fundamental origin of some widely applicable scaling laws, as the Gutenberg-Richter frequency-magnitude relationship, the modified Omori law or the Omori-Utsu law for aftershock decay and Bath's law for the difference between the magnitude of the largest aftershock and mainshock.
The frequency-magnitude distribution (Gutenberg and Richter, 1944), firstly examinated, describes the relation between the frequency of occurence and magnitude of earthquake 10 ( ) log N m a bm = - (1) where, N (m) is the cumulative number of events with magnitude larger than or equal to m a and b are constants. The parameter a showing the activity level of seismicity exhibits significant variation from region to region becauses it depends on the period of observations and area of investigation. The parameter b describes the size distribution of events and is related to tectonic characteristics of the region under investigation. It is shown that the estimated coefficient b varies mostly from 0.6 to 1.4 (Weimer and Katsumata, 1999). Many factors can cause perturbations of the b value. For example it is established that the least square procedure introduce systematic biais in the estimation of the b parameter (Marzocchi and Sandri, 2003;Sandri and Marzocchi, 2005). The b value of a region not only reflects the relative proportion of the number of large and small earthquakes in the region, but is also related to the stress condition over the region. The physical implication of the b value, however, is not as obvious, is still under investigation.
The temporal decay of aftershocks is important task because it contains information about seismogenic process and physical condition in the source region. It is well known that the occurrence rate of aftershock sequences in time is empirically well described by the Omori-Utsu law as proposed by Utsu (Utsu, 1969). The power law decay represented by the Omori-Utsu relation is an example of temporal self-similarity of the earthquake source process. The variability of the parameter p-value is related to the structural heterogeneity, stress and temperature in the crust (Mogi, 1962;Kisslinger and Jones, 1991). The importance to derive reliable p-value, is due to the fact that it contains information about the mechanism of stress relaxation and friction law in seismogenic zones (Mikumo and Miyatake, 1979), but this information cannot be derived without a precise characterization of the empirical relations that best fit the data. In this study we examine the temporal patterns of aftershock distribution of earthquakes occurred in Algeria-Morocco region. The parameters of Omori-Utsu law are estimated by the maximum likelihood method, assuming that the seismicity follow a nonstationary Poisson process. Two statistical models, a first stage Omori-Utsu model, model without secondary aftershock and a two stage Omori-Utsu model including the existence of secondary aftershock activity, are tested for goodness of fit to aftershock data. We adopt the Akaike Information Criterion denoted AIC, (Akaike, 1974) as a measure for selecting the best among competing models, using fixed data.
A model describing the temporal behavior of cumulative moment release in aftershock sequences is analyzed as alternative approach with respect to the Omori-Utsu law. In this model the static fatigue is assumed to be the principal explanation of the aftershock temporal behavior. Following Marcellini (1997), the different aftershock sequences are analyzed using this model. The most important result shows that globally the model explain and fit correctly the data.
Modified Bath's law is analyzed for each aftershock sequence, it allows us to derive some important properties related to the energy partionning. This approach is used to derive the energy dropped during the mainshock.
The fractal dimension D 2 deduced from the correlation integral is used to carry out the spatial analysis of the studied sequences. The inter-event distances analysis is performed for the different aftershock sequences and the density of probability of the inter-event distances are derived using non-parametric approach, especially using the kernel density estimation technics (Silverma, 1986).

Tectonic sketch and seismicity overview
In this section we give a large overview of the tectonic skech and the seismicity of the studied region (   , moreover reverse fault, many of which has N70ºE to E-W direction. In many cases, the cited strike-slip faults moved mainly as normal faults, releasing by this way the regional tension, practically perpendicular to the compression. In this study, we analyze the aftershock sequences triggered by some important and damaging earthquakes which occurred in the Morocco-Algeria region. The seismic activity in this region is mainly characterized by moderate to destructive magnitude events. It has been the site of numereous historical earthquakes, and was therefore subject to extensive damage and several thousands of casualties in the past. It is well known that the process to identify aftershock is related to the seismicity declustering process, a crucial step in separating an erthquake catalog into foreshock, aftershock and mainshock. This process is widely used in seismology, in particular for seismic hazard assessment and in earthquake prediction model. There are sevceral declustering algorithms that have been proposed. Up to now, most users have applied either the algorithm of Gardner and Knopoff (1974) or Reasenberg (1985). Gardner and Knopoff (1974) introduced a procedure for identifying aftershocks within seismicity catalog using interevent distances in time and space. They also provided specific space-time distances as a function of mainshock magnitude to identify aftershocks. This method is known as a window method and is one of the most used. Reasenberg's algorithm (1985) allows tolink up aftershock triggering within an earthquake cluster: if A is the mainshock of B, and B is the mainshock of C, then all A, B and C are considered to belong to one common cluster. When defining a cluster, only the largest earthquake is finaly kept to be the cluster's mainshock. In this study we have used the single-link cluster algorithm (Frohlich and Davis, 1990). We first identify large magnitude within the database. Earthquakes occurring within a certain distance d in km and periode T in days after these potentialmainshock were extrated from the database. Aftershocks were then selected from this subset of events using single-link cluster algorithm (Frohlich and Davis, 1990). A space-time metric D is used to define the proximity of events relative to one another. The metric is defined as where d is the distance between two epicenters Δ t is the temporal separation between two events and B is a constant relating time to distance. Afteshpocks are then defined as events within a D c -size (critical size) cluster containing the mainshock. The single-link cluster analysis has been used with the parameters defined by Frohlich and Davis (1990), it is natural choice for aftershock selection because it avoids subjective decision-making with regard to the spatial distribution and duration of the sequence, it requires no assumption as the efficiency of event triggering, as is necessary when applying other clustering algorithms (e.g. Reasenberg, 1985).
The May 26, 1994 Al Hoceima earthquake (Mw 6.0) took place at 12 km depth, the focal mechanism indicates the presence of a main set of sinistral fault with a N-S trend, which may involve several parallel surface (Calvert et al., 1997 the region is affected by several faults. The regional seismicity analysis shows that the Babors chain seems to belongs to a "transition zone" between a large belt of reverse faulting along the western and central part of northern Algeria and a more distributed zone where deformation is mainly accomodated through strike-slip faulting (Beldjoudi et al., 2009). The aftershock sequence include 111 of the best recorded events with more than 54 with RMS < 0.3s, ERH < 3km and ERZ < 3km. The event magnitude varies between 1.3 and 5.2.

Magnitude-frequency relationship
The frequency-size distribution (Gutenberg and Richter, 1944) describes the relation between the frequency of occurrence and the magnitude of earthquake Many factors can bias the b value estimation. It is established for example that the least square procedure introduce systematic error in the estimation of the b parameters (Marzocchi and Sandri, 2003;Sandri and Marzocchi, 2005). Locally, the b value not only reflects the relative proportion of the number of large and small earthquakes in the region, but it is also related to the stress condition over the region (Mogi, 1962). The physical meaning of the b value is, however not clear and still under investigation.
The first step in the analysis of the Gutenberg-Richter law, is the determination of the threshold magnitude of completness. It is defined as the lowest magnitude for which all the events are reliably detected (Rydelek and Sacks, 1989). There are many approaches to the estimation of m c , the most popular are; The maximum curvature (MAXC) method (Weimer and Wyss, 2000) defines the completeness magnitude as the magnitude for which the first derivative of the frequency magnitude curve is maximum (being the maximum of the non-cumulative frequency-magnitude distribution). The goodness of fit (GFT) method (Weimer and Wyss, 2000) compares the observed frequencymagnitude distribution with a synthetic distribution, and the goodness of fit is calculated as the absolute difference of the number of earthquake in each magnitude bins between the observed and synthetic distribution. The synthetic distribution is calculated using a and bvalues estimated from the observed dataset by increasing the cutoff magnitude. The completness magnitude, is given by the magnitude for which 90% of the data are fitted by a straight line. The entire magnitude range (EMR) method was developed by Ogata and Katsura (1993) and modified later by Woessner and Weimer (2005). The maximum likelihood estimation method is used to estimate the power law G-R law parameters a and b. The same method is applied to estimate the mean and standard deviation of the Normal distribution considered for the incomplete part of the frequency-magnitude distribution. μ and σ are the magnitude at which 50% of the earthquakes are detected and the standard deviation describing the width of the range where earthquakes are partially detected, respectively.
The frequency-magnitude distribution as shown previously is defined as, (Gutenberg and Richter, 1954) where, N ( ≥ m) is the cumlulative number of earthquakes with magnitudes greater than or equal to magnitude m, and a and b are constants. This relationship holds for global earthquake catalogs and is applicable to aftershock sequences as well. The estimation of the b-value has been the subject of considerable research, and various methods exist. The most statistically appealing of these is the maximum likelihood method first used independently by Aki (1965) and Utsu (1965), which gives the estimate of the b-value as where m min is the minimum magnitude up to which the Gutenberg-Richter law can accurately represent the cumuluative number of earthquakes larger than or equal to a given magnotude, and m is the average magnitude. The definition of m min origionates from the definition of the probability density function for magnitudes, f (m), consistent with the mathematical form of the Gutenberg-Richter law. This is given as (Bender, 1983) where m max is the maximum magnitude up to which f (m) describes the distribution of magnitudes and λ = bLn (10). There are two problems with equation (4) when applied to real earthquake catalogs. Firstly, it considers reported earthquake magnitudes as a continuous variable, which is not accurate as most earthquake catalogs report magnitudes up to a precision of one decimal place for each event. Therefore, magnitude should be considered as a grouped or binned variable with a finite non zero bin length δm, which generally and in our case it is taken equal to 0.1, δm = 0.1. This means that instead of observing m min physically in the catalog, we observe a different minimum magnitude m c = m min + δm 2 , which is the aforementioned completeness magnitude. This is the smallest observed magnitude in the catalog above which the cumulative number of earthquakes, above a given magnitude, are accurately described by the Gutenberg-Richter law. It is very important to note that the minimum magnitude of completeness m c , physically observed in the catalog, is not present in the equation (4), and hence we need to express m min in terms of m c to be able to use the formula for real catalog.
Assuming that m min = m c would lead to a bias in the estimate. The other source of inaccuray in equation (4) is the fact that it considers the maximum magnitude value for the dataset to be infinite, which is never the case. in fact, for aftershock sequences, the maximum magnitude in the sequence can be pretty small. This too introduces a bias in the estimate of the b-value obtained using the equation (4) (Bender, 1983;Tinti and Mulargia, 1987;Guttorp and Hopkins, 1986). This problem was considered by Bender (1983), and she gave a method for obtaining the b-value of grouped and finite maximum magnitude data as a function of the root of the following equation where q = exp − Ln (10).b .δm and n = (m max − m min ) δm .
It is straightforward to see that for n → ∞, equation (6) gives (Tinti et al.,1987;Gutorp et al., 1986). 10 1 (Tinti and Mulargia, 1987;Gutorp and Hopkins, 1986). It has been observed, however, that for a difference of about 3.0 in magnitude between m min and m max , the value of b obtained from equation (6) agrees closely with the asymptotic value of b obtained from equation (7) (Bender, 1983). In Fig.2  Laalam 2006 sequence. The higher value obtained for El Asnam series reflect the quality of the data include in the El Asnam series, it seems that the file doesn't include all the events. It is till now one of the "best" file including the mainshock and almost all the aftershock occurred just after the mainshock. The b-value obtained for these two series, 0.82 ± 0.10 for El Asnam series and 0.99 ± 0.10 for Laalam aftershock series are close to 1.0, a typical value for the aftershock sequences. For the Zemouri aftershock sequence we have used the best combination between m c obtained for 95% and 90% confidence and the maximum curvature procedure (Weimer and Wyss, 2000), which gives a threshold magnitude equal to m c = 3.5 m bLg . To improve the b-value analysis for this sequence, we obtained a threshold magnitude m c = 3.7 m bLg , using the procedure based on the stabilisation of the b-value and its uncertainty derived by Shi and Bolt (1982), as implemented by Weimer and Zuniga, (1994) in the software Zmap under Matlab. We obtained b = 1.10 ± 0.04 for the first procedure and b = 1.30 ± 0.06 for the second procedure, the difference between the two values is of the order of 0.20.

Decay rate of aftershock activity -Omori-Utsu law
The third studied scaling law is related to the decay rate of aftershock activity. It is well know that the decay rate of aftershock activity with time is governed by the modified Omori law or Omori-Utsu law (Utsu et al. 1995), where, θ = (k , p , c) are the model parameters.
The log likelihood is then given by;  In our case, Fisher information matrix is given by the results obtained using this procedure for the aftershock sequences are shown on Fig. 4.
It is often observed that a sequence of aftershocks contains secondary aftershocks, which are aftershocks of a major aftershock (Utsu, 1970). Secondary aftershock are typically detected as changing-point in the activity rate of the sequence using Akaike information criteria (AIC) (Akaike, 1974). In this study, the changing-point in the activity rate of the sequence is detected by using the plot of cumulative number of events vs time from the mainshock. Assuming one secondary aftershock occurred at time τ 0 , we test four different hypothesis, H1 : no secondary aftershock, H2: a secondary aftershock series does exist with the same p-value and c-value as the original series, H3; secondary aftershock series does exist with same c-value as the original series and H4: a secondary aftershock series does exist with Omori law parameters different from the original series. Four point process models are tested. Respectively, the conditional intensity and the cumulative function of such point process are given by; and the cumulative function is given by; T s = 0 is the occurrence time of the main shock, and T e is the occurrence time of the last event.
All the occurrence time are counted as the number of days elapsed from the mainshock. For all the model relative to the hypothesis H1, H2, H3 and H4, the model with the lowest AIC is selected as the best model. The AIC is used as a measure to select which model fits the observations better. This is a measure of which model most frequently reproduces features similar to the given observations, and is defined by  Table 1.
As pointed previously, it is often observed that a sequence of aftershocks contains secondary aftershocks, aftershocks of a major aftershock (Utsu, 1970). If a secondary aftershock equence starts at time τ 0 then from the Eq. 15 and 16, the conditional intensity of the point process is given as    from the graphs of the residual process we can deduce, as shown previously using the AIC criteria, a global agreement with the Omori-Utsu model. We proceded with the same technics for the others aftershock sequences, the AIC is used to select the most appropriate model fitting the data. Table 2, gives the obtained results.   The temporal aftershock decay is also analyzed using the approach introduced by Marcellini (1997). This approach describes the temporal behavior of the cumulative seismic moment released in aftershock sequences. It is an alternative approach to the Omori-Utsu law, previously analyzed. Static fatigue is assumed to be the principal explanation of the aftershock temporal behavior. Under the condition that the main shock causes a redistribution of stress, the initial stress condition of the afterhock sequence at main shock origin time t 0 can be considered as the superposition of the stress before the main shock and the stress step dσ caused by the dynamic rupture of the main shock. It has been shown Marcellini (1995) that ; ; V k is the cumulative focal volume, V 0 is the focal volume of the main shock, V j is the focal volume of the j th aftershock and M 0m and M 0 j are the seismic moment of the main shock and the j th aftershock, respectively.
Following Marcellini (1997), we consider the aftershock zone as characterized by barriers that breaks after a given elapsed time proportional to the stress intensity factor K i and therefore to the stress σ i . Data fit of equation (22) characterizes the distribution, namely, the more our data explained by Eq. 22, the most likely the distribution of σ is close to the uniform distribution. This property is partially influenced by the spatial variation of the stress change produced by the mainshock, given that in the present static fatigue model a barrier breaks at σ i , which is the superposition of the static stress before the mainshock and stress step induced by the mainshock In this study, as in Marcellini (1997), the seismic moment is evaluated using the Thatcher and Hanks (1973) relation, given by to test if the static model fatigue holds as represented by equation (22), two conditions must be checked:  Taking into account the obtained coefficient of correlation r 2 and given on each graph, the modele y = a + b Log(t) introduced by Marcellini (1997) fit very well our data.

Energy partitioning
The other scaling law exalined in this study is the modified Bath law. In its original form, Bath law states that the difference Δm between a given mainshock with magnitude m ms and its largest aftershock magnitude m as max is approximately constant, independently of the mainshock magnitude (Bath, 1965). That is, we derive the partionning of the energy, in the following way. It is wellestablished that the magnitude distribution of aftershocks clearly exhibits a near-universal scaling relative to the mainshock magnitude. To explore this relation for our aftershock sequences, we will determine the ratio of the total energy radiated by the afytershock sequence to the seismic energy radiated by the mainshock. The energy radiated during an earthquake is related empirically to its moment magnitude m by (Utsu, 2002) the ratio of the total radiated energy by the aftershocks E as to the total energy radiated by the mainshock E ms is given by ( )  assuming that all earthquakes have the same seismic efficiency, which means that the ratio of the radiated energy to the total drop is stored as elastic energy is also the ratio of the drop in the stored elastic energy due to the aftershocks to the drop in the stored elastic energy due to the mainshock. Finally, the following relation is derived Using the Eq. 39 for the studied aftershock sequences, the results are shown on the following  Table 3. b-value, Bath law and energy partitionning derived for the different studied aftershock sequences.
From the obtained results shown on Table 3, we deduce the percentage of the total energy radiated during the mainshock. From this point of view, 95 % of the total energy has been radiated during the Al Hoceima 1994 mainshock and 65% during Al Hoceima 2004 mainshock.
On the other side the El Asnam 1980 main shock radiated about 98% of the total energy, 2% has been radiated by the aftershock sequence, but we shouldpoint out that these results depend directly on the quality of data used and it is clear that whatever the sequence of aftershocks used, it is still incomplete, especially at the beginning just after the occurrence of the main shock. For the aftershock sequence triggered by the El Asnam earthquake of 1980, it seems that the sequence used is truncated due to the delay in the implementation of seismological network just after the main shock. Nevertheles, the results shown on Table 3, gives a large overview on the ratio of the total energy radiated by the main shock and by the aftershocks. Thus, 86 % of the total energy has been radiated during the main shock of Zemouri 2003 and 99% during Laalam 2006 mainshock.

Spatial aftershock distribution
It is well known that seismicity is a classical example of a complex phenomenon that can be quantified using fractal theory (Turcotte, 1997). In particular, fault networks and epicenter distributions have fractal properties (Goltz, 1998). Thus, a natural way to analyze the spatial distribution of seismicity is to determine the fractal dimension D 2 . This D 2 -value is an extension of the Euclidian dimension and measures the degree of clustering of earthquakes.
In the two dimensional space, D 2 can be a decimal number and ranges from 0 to 2.0. Therefore In practice, however, for large values of r the gradient is artificially low, whereas for small values of r the gradient is artificially high. These two cases have been called "saturation" and "depopulation" (Nerenberg and Essex, 1990). Whereas, it is common in the estimation of the fractal dimension to use fiting procedure to straight line than to a subjectivelly chosen straight part of the curve. Nerenberg and Essex (1990), give formulas for determining the distances of depopulation and saturation, r d and r s given by where d is the dimensionality of the data cluster, and 2R is the approximate lenght of the side of the area containing data. As pointed by Eneva (1996), it is safe to start the scaling range at values of r as low as r d / 3, but in this study, we have use the slope method, which as explained previously consist in estimating the slope of the double logarithmic plot of the correlation integral versus distance. The stability of the scaling range is verified with the slope. It consists in calculating the first derivative between each two pints of the correlation integrale curve and plotting versus logarithmic distance. To eliminate the depopulation and saturation effects, the scaling range is defined within the region where the slope is most constant. Using this procedure and as for Al Hoceima 1994 and 2004 series, the results obtained are displayed in the following figure The results obtained are close to 2.0, which allow us to deduce that the spatial distribution of the epicenters tends to be uniform on the plane.
The ratio of the slip on the primary fault to the total slip over the fault system is given by (Khattri, 1995) ( ) where S p is the slip over primary fault and S represents the total slip over the fault-system.
The obtained results are shown on  Table 4. Fractale dimension D 2 and ratio of the slip on the primary fault over the fault system.
In this section we attempt to analysis the inter-event distance distribution of probability. We use a non-parametric approach to analysis the density of probability of the inter-event distances, especially the kernel density estimation, this methodlogy is clearly presented in Silverman (1986). We used it in the following way. Given a sample of n observations x 1 , x 2 ,........., x n with unknown probability distribution function f , the kernel density estimate of f is given by, where K (.) is a positive function called kernel, a typical example is the Gaussian kernel defined by; ( ) in Eq. 44, h is a parameter and its determination is crucial. There are differents method to estimate the parameter h , in the following lines we summarise the least-square cross-validation introduced by Silverman (1986, pp 48). The kernel density estimate of f is given as shown previously by considere a sample of n observations x 1 , x 2 ,........., x n ; the quadratic uncertainty is then given The last term of the right member of the last equality is independnat of the parameter h , thus the optimal value of h is obtained by minimizing the two others terms. The problem then consists to find h minimizind the score function defined as follow futhermore, the score function could be written in the following form with, (51) assuming that the minimum of M 0 (h ) is in the vicinity of the minimum of E M 0 (h ) and then Also, for n large, (n − 1) ≅ n; finaly the problem consites to minimize the simplified score function defined as follow, we replace M 0 (h ) by (52) we observe that in the case of the Gaussian kernel, the kernel K (2) (.) is a centered Gaussian distribution with variance 2. In this case the simplified score function is written as follow the parameter h minimizing the simplified score function is then obtained as solution of the following equation the parameter h obtained using teh cross-validation method is then given by where Proof. The proof of this theorem is given in details in Silverman (1981).
Under the null hypothesis, samples can be simulated from the kernel density estimate by using  is a randomly simulated sample from the standard normal distribution.
if h crit * (k) is the parameter obtained from x _ then the p value of this test is given by # is the sign of the number of element in the set. In practice, we apply the series of tests  The Fig. 11(a) gives the estimated density obtained using the rule of thumb. The cross validation optimal smoothing method gives a bandwidth parameter h = 1.2, the estimated density of probability obtained is shown on Fig. 11 The estimated density using the cross validation optimal smoothing has been obtained with a bandwidth parameter h = 1.0. The density displays an inter-event distances concentration around the value 4.29 km as shown on Fig. 13(b). The estimated density with one mode critical bandwidth h crit = 6.6; displays a concentration of the inter-event distances around 9.0 km.

Conclusion
Aftershock sequences in Algeria-Morocco region have been analyzed in order to estimate and derive with accuracy the parameters of the most important scaling laws in statistical seismology. For each aftershock sequence the threshold complteness magnitude m c has been derived using different approach, we have choose the one giving the "most" stability of the fit of the cumulative number by a straight line. For the magnitude above the threshold magnitude above the threshold completeness m c , we used the maximum likelihood to estimate above the threshold completeness m c . The results obtained are close to 1.0, the typical universal value for aftershock sequence. The Omori-Utsu law for aftershock decay and Bath's law for the difference between magnitude and the largest aftershock and mainshock as modified by Shcherbakov et al., (2004a, b)., while using AIC as a measure to select the most appropriate model between the first stage and two stage Omori-Utsu model, the rate of decay of aftershocks was found to follow in all case, except the May 21, 2003 Zemouri aftershock sequence, the first stage Omori-Utsu law, which mean the Omori-Utsu model without secondary aftershock. The decay of aftershock trigged by the May 21, 2003 Zemouri earthquake exhibit a better fit with a two stage Omori-Utsu model, denoted model 4 in the text. The later allow us to include in the modelisation the secondary aftershock of magnitude 5.8 m bLg . The study of the Omori-Utsu law has been complemented of aftershock sequences. This procedure is an alternative to the Omotri-Utsu law based on the physics of the static fatigue mechanism. The study of Gutenberg-Richter relationship and the modified Bath's law, allowed us to perform the partitioning of the energy released during each aftershock sequence by the main shock and aftershocks. The spatial analysis has been performed using the fractal dimension D 2 derived using the Correlation integrale. Using non parametric estimation approach, especially the kernel density estimation, we derive the estimation of the density of the probability of the inter-event distances. The mode of this density highlight how the spatial clusters of the events.