Open access

Scaling Properties of Aftershock Sequences in Algeria-Morocco Region

Written By

M. Hamdache, J.A. Pelàez and A. Talbi

Submitted: 04 May 2012 Published: 20 March 2013

DOI: 10.5772/54888

From the Edited Volume

Earthquake Research and Analysis - New Advances in Seismology

Edited by Sebastiano D'Amico

Chapter metrics overview

2,175 Chapter Downloads

View Full Metrics

1. 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 (Mw6.0) and February 24, 2004 (Mw6.1) which occurred in northern Morocco, the October 10, 1980 El Asnam earthquake (Mw7.3), the May 21, 2003 Zemouri earthquake (Mw6.9) and March 26 2006 Laalam earthquake (Mw5.2) in northern Algeria.

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


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 non-stationary 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 D2 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).


2. Tectonic sketch and seismicity overview

In this section we give a large overview of the tectonic skech and the seismicity of the studied region (Pelaéz et al., 2012). The Morocco-Algeria region, namely the Maghrebian region (Fig. 1) occupies the NW part of the African (Nubia) Plate in what is referred to its continental crust. Its oceanic crust continues till the area of the Azores Islands. To the North it is immediately situated the Eurasian Plate, although between the Gibraltar Arc and the South of Italy an intermediate complex domain is intercalated. This domain is formed by some oceanic basins, as is the Algero-Provençal Basin and the Thyrrehnian Basin, and by a former region, presently disintegrated and now forming the Betic-Rifean Internal Zone, the Kabylias (in Algeria), the Peloritani Mountains (Sicily) and the Calabrian area in Italy (AlKaPeCa domain). This area underwent from the early Miocene a northwards subduction of Africa, then opening the small oceanic basins quoted, accompanied by the disintegration of the AlKaPeCa domain.

Figure 1.

Tectonic sketch showing the main tectonic domains in the studied region.

Presently, the convergence between the Nubia Plate and Iberia has an approximate NNW-SSE direction, with values of the order of 3 to 5 cm/year, according the places. This compression is accompanied, at least in the area of the Gibraltar Arc (in the Alboran Sea) by a noticeable ENE-WSW tension, in some cases even more important than the compression. For this reason, in the Alboran area, some extensional movements can be important The Maghrebian region is also a complex area in which the Saharan Shield affected by the Pan-African Orogeny (Precambrian to early Cambrian) is in contact with the Atlasic Mountains of mainly Alpine age (Fig. 1). To the North of the Atlas is situated the Moroccan Meseta, to the West, and the High Plateaus in Algeria, which to the North contact with the Rif and Tell mountains, typically Alpine chains. The Saharan Shield forms part of the Precambrian areas of Africa, clearly cratonized and generally not affected by later important deformations. In fact, in the Maghrebian area it corresponds to a clearly stable area. In Morocco, the so called Antiatlas, corresponds to a Precambrian and, mainly, Paleozoic area, making a tectonic transition between the shield and the Atlas. The Atlasic Mountains correspond to an intracontinental chain. To the West, in Morocco, the High Atlas reach the coast in the Agadir area and continues to the Northeast and East, passing, although with lesser heights, to the Saharan Atlas, which cross Algeria and reach the central part of Tunis. They can be considered as aulacogens bordering the northern part of the Saharan shield. To the North, the Middle Atlas in Morocco has a different direction, NE-SW, separating the Moroccan Meseta and the High Plateaus, both forming by Paleozoic rocks, although with a Mesozoic and Tertiary cover, well developed in some areas. On the whole, the Atlasic Mountains has been tectonically unstable from the Triassic times, and along the Alpine orogeny suffered important deformations and, more recently, also important volcanism, reaching the Quaternary. The Rif and Tell thrust southwards the Moroccan Meseta and the High Plateaus, and even in some places part of the Atlasic Mountains. They are formed by sedimentary External zones (only slightly affected by metamorphism in some Moroccan places) and by Internal zones. Mostly Internal zones (divided in several tectonic complexes) are affected by alpine metamorphism, moreover the existence of previous Paleozoic and even older deformations. In any case, their present structure has being formed during the Alpine Orogeny. They appear mainly to the E of Tetuan, in Morocco, and in the Kabylias, in Algeria. These Alpine chains have being structured from the Cretaceous to the Oligocene-early Miocene. Later, were formed numerous Neogene basins, clearly cutting in many cases previous structures. In this time, particularly from the late Miocene to the present, a near N-S compression provoked the existence of strike-slip faults (NE-SW, sinistral, and NW-SE, dextral), 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. Among the largest recorded earthquake in Morocco, the February 29, 1960 Agadir (Morocco) earthquake (Ms6.0), the May 26, 1994 and February 24, 2004 Al Hoceima earthquakes with Mw6.0 and Mw6.1 respectively (Pelàez et al., 2007). In northern Algeria, earthquakes up to magnitude Ms 7.3 have been recorded, namely the October 10, 1980 El Asnam earthquake (Ms 7.3) (Ouyed et al.,1981) and May 21, 2003 Zemouri earthquake (Mw6.9) earthquake, (Hamdache et al., 2004). It is important to point out that these two events were the most damaging events occurred in northern Algeria and, the region of EL Asnam experienced in the past damaging earthquake on 09 September 1954 (Hamdache et al., 2010). The original epicentral database and the arrival time readings of the different aftershock sequences used, were obtained from the Spanish National Geographical Institut (IGN) earthquake catalog. For the sequences triggered by the May 21, 2003 Zemouri earthquake and by March 20, 2006 Laalam event, we used the data recorded by portable seismological stations network monitored by CRAAG. 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 inter-event 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 Dc- 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). Analysis of the aftershock sequence highlighted the presence of NNE-SSW distribution of the seismicity (El Alami et al., 1998). The sequence used covers a period of about one year from the mainshock, including 318 events with magnitude ranged from 2.0 to 6.0. The February 24, 2004 Al Hoceima earthquake (Mw 6.1) took place at a depth of 10 to 14 km. The focal mechanism suggests that the active nodal plane corresponds to a sinistral strike-slip fault oriented N11 N and dipping 70 toward the E (Stich et al., 2005). The aftershocks were aligned preferently in NNE-SSW, in the same way as one of the nodal planes of the focal mechanism. As pointed by Galindo-Zaldivar et al., (2009) the main stresses for the two aftershock sequences trigged by the May 26, 1994 and February 24, 2004 events have a trend with NW-SE compression (P-axis) and orthogonal NE-SW extension (T-axis) compatible with the present convergence of the Africa and Eurasia plate. The aftershock sequence includes 1233 events with magnitude ranged between 0.6 to 6.1. The time span of the sequence is about 1 year.

The October 10, 1980 El Asnam earthquake (Mw 7.3) is one of the most important and most damaging event occurred in northern Algeria. This event has taken place on the Oued-Fodda reversse fault. The later is segmented into three segments, ruptured along 26 km. This fault is located in the Cheliff high seismogenic Quaternary bassin, considered as very active. It is important to point out that almost all the seismicity in northern Algeria is located around the Plio-Quaternary intermountains active bassins (Meghraoui et al., 1996). The aftershock sequence used include the 130 most important magnitude events, ranged between 2.4 and 7.3. In the same way, the May 21, 2003 earthquake (Mw 6.9) has been located in the NE continuation of the south reversse fault system of the Quaternary Mitidja bassin (Maouche et al., 2011). The aftershock sequence we used include 1555 magnitude events, ranged between 0.9 to 6.9 and recorded during the first 40 days from the main shock (Hamdache et al., 2004). The March 20, 2006 Laalam earthquake (Mw 5.2), was located in the Babor chain, in the 'Petite Kabylie' south of Bejaia city. This chain belongs to the Tell Atlas, which is a portion of the Alpine belt in northern Africa. As pointed out by Beldjoudi et al., (2009) 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.


3. Magnitude-frequency relationship

The frequency-size distribution (Gutenberg and Richter, 1944) describes the relation between the frequency of occurrence and the magnitude of earthquake


the statistic N(m) is the cumulative number of events with magnitude larger than or equal to m, whereas a and b are positive constant. The parameter a which exhibits significant variations in space measures the activity level of seismicity. This parameter is sensitive to the input period of observation and area of investigation. The parameters 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 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 mc, 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 frequency-magnitude 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 b-values 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 mmin 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 mmin 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 mmax 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 mmin physically in the catalog, we observe a different minimum magnitude mc=mmin+δm2, 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 mc, physically observed in the catalog, is not present in the equation (4), and hence we need to express mmin in terms of mc to be able to use the formula for real catalog. Assuming that mmin=mc 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=(mmaxmmin)δm.

It is straightforward to see that for n, equation (6) gives (Tinti et al.,1987; Gutorp et al., 1986).


(Tinti and Mulargia, 1987; Gutorp and Hopkins, 1986). It has been observed, however, that for a difference of about 3.0 in magnitude between mmin and mmax, the value of b^ obtained from equation (6) agrees closely with the asymptotic value of b^ obtained from equation (7) (Bender, 1983).

For the two aftershock sequences triggered by the mainshocks occurred around Al Hoceima city in Morocco on 1994 and 2004, the threshold magnitude mc has been examined in details, using the different procedures introduced previously, especially the maximum curvature procedure MAXC (Weimer and Wyss, 2000) and the changing point procedure introduced by Amores, (2007). The results obtained for these two aftershock sequences, are shown in Fig. 2

Figure 2.

Graphs showing the non cumulative number of events, the threshold magnitide and the adjustment of the cumulative number by straight line with equation Log10N(m)=abm for mMc. (a) For Al Hoceima 1994 afeshock sequence and (b) for Al Hoceima 2004 aftershock sequence.

In Fig.2, the frequency-magnitude relation for the 1994 and 2004 aftershocks series of Al Hoceima (Morocco) are displayed. Based on maximum curvature procedure (MAXC), the magnitude of completeness was taken equal to 2.8 for the 1994 aftershock seqiuence and 3.4 for the 2004 sequence. It is important to point out that the changing point procedure (Amores, 2007) gives the same results. Using these threshold magnitudes, we derive the b-value of the Gutenberg-Richter relationship and its standard deviation using the maximum likelihood procedure. The b-value is estimated equal to 0.92±0.02 for the 1994 aftershock sequence and 1.073±0.003 for the 2004 series. For both sequences the b-value obtained is close to 1.0, the typical value for aftershock sequence. The obtained parameter a is equal to 4.689±0.058 for the 1994 aftershock series and equal to 6.305±0.014 for the 2004 aftershock sequence. The results obtained for the aftershock sequences triggered by the events occurred in northern Algeria, namely the October 10, 1980 El Asnam earthquake (Mw 7.3), the May 21, 2003 Zemouri earthquake (Mw 6.9) and March 20, 2006 Laalam earthquake (Mw 5.2) are displayed in the following Fig. 3.

Figure 3.

Graph displaying the non cumulative number of event in blue circle and the adjustment of the cumulative number (red circles) by straight line, representing the Gutenberg-Richter relation for the aftershock sequence triggered by (a) October 10, 1980 El Asnam earthquake (b) by May 21, 2003 Zemouri earthquake and (c) by March 20, 2006 Laalam earthquake.

For the sequence of October 10,1980 El Asnam earthquake and March 20, 2006 Laalam earthque, we have used the maximum curvature procedure to derive the threshold magnitude, we obtained mc=3.9mbLg for the El Asnam 1980 aftershock sequence and mc=2.1mbLg for 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 mc obtained for 95% and 90% confidence and the maximum curvature procedure (Weimer and Wyss, 2000), which gives a threshold magnitude equal to mc=3.5mbLg. To improve the b-value analysis for this sequence, we obtained a threshold magnitude mc=3.7mbLg, 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.


4. 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),


λ(t) is the rate of occurrence of aftershocks per unit time, at time t after the mainshock (t=0). The parameters k,c and p are positive constants. k depends on the total number of events in the sequence, c on the rate of activity in the earliest part of the sequence and p is related to the power law decay of aftershocks. It is generally accepted that the number of aftershocks cannot be counted completely in the beginning of the sequence when small shock are often obscured by large ones due to overlaping, thus an overly large value of the parameter c is obtained. After Utsu, (1971), the parameter c must be zero if all shocks could be counted. Two opinions constitue the nowdays debate around the c value: one is that c value is essentially zero and all reported positive c value results from incompletness in the early stage of an aftershock sequence; the second point of view is that c value do exist (Enescu and Ito, 2002; Kagan, 2004; Shcherbakov et al. 2004a, b; 2005, 2006; Enescu et al. 2009). The constant c is a controversial quantity (Utsu et al. 1995; Enescu et al. 2009) and is mainly influenced by the incomplete detection of small aftershocks in the early stage of the sequence (Kisslinger and Jones, 1991). According to Olsson (1999), p values generally vary in the interval 0.5 to 1.8 and this index has shown by Utsu et al. (1995), differs from sequence to sequence and vary according to the tectonic condition of the region, however it is not clear which factor is the most significant in controlling p value. Although more attention in the estimation of the p value have been done to take into account the recommandation by Nyffenegger and Frolich (1998, 2000), we have use the standard use the standard deviation obtained by maximum likelihood which could be over or underestimating (Nyffenegger and Frolich, 1998; 2000).

In this study aftershock sequences are modeled using point process defined by the following conditional intensity,


with, t is the history of observation until time t, e.g the σ- algebra generated by the events occurred before the time t, thus, following Daley and Vere-Jones (2005), the family (t,t+) is called the natural filtration of the point process or the internal history. The conditional intensity in equation 4, is independent of the internal history and depends only on the current time t, like λ(t|t). It defines a non-stationary (nonhomogeneous) Poisson process. We assume the occurrence time's of the aftershock sequence, namely {t1,t2,.............,tn}[Ts,Te] are distributed according to a non-stationary Poisson process with conditional intensity function defined by equation 4. The parameters k,p and c are estimated by maximizing the following likelihood function


where, θ=(k,p,c) are the model parameters.

The log likelihood is then given by;




it follows that the maximum likelihood estimate MLE θ˜=(k˜,p˜,c˜) of the parameter θ is solution of the following equation


Following, Ogata (1983), the maximum estimation of the parameters are obtained by using the Davidon-Fletcher-Powell optimization algorithm (Press et al. 1986, pages 277, 308) applied to equation 14. The standard deviation of the estimated parameters through the maximum likelihood procedure are derived using the inverse of the Fisher information matrix J(θ˜)1. 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;


Ts=0 is the occurrence time of the main shock, and Te 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


a model with smaller value of AIC is considered to be better fit to the observations. The fit of the data by the Omori-Utsu law is obtained for both the minimlum magnitude mmin in the aftershock sequence and the threshold magnitude mc discussed and derived in the previous section. The results of the fit are shown on the following figure

Figure 4.

Fit of the number of events pêr day by the modified Omori law, for the studied aftershock sequence. The graphs display the fit for mmc in dashed line and in solid line for mmmin.

The Omori-Utsu parameters obtained for each aftershock sequence for magnitude above the threshold magnitude, mmc, are given on the Table 1.

Omori-Utsu parameters for mmc
p±σp c±σc k±σk
Al Hoceima 1994 0.76 ± 0.03 0.0040 ± 0.0008 12.47 ±1.51
Al Hoceima 2004 0.85 ± 0.02 0.034 ±0.016 47.38 ±3.81
El Asnam 1980 0.84 ±0.06 0.025 ±0.045 4.48 ±3.81
Zemouri 2003 1.13 ±0.06 0.311 ±0.053 107.53 ±15.83
Laalam 2006 0.69 ±0.09 0.014 10.27 ±1.47

Table 1.

Omori-Utsu Parameters p,c,k and their standard deviation obtained for each aftershock sequence with magnitude above the threshold magnitude, mmc

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


denoting, θ=(k1,p1,c1,k2,p2,c2,τ0) and where H(t) is teh Heaviside unit step function. The occurrence of secondary aftershock is analyzed using the procedure of detecting the aftershock activity change point by AIC, procedure introduced by Ogata (1999), Ogata (2001) and Ogata et al., (2003). Using the changing point AIC procedure to the aftershock sequence of Al Hoceima 1994, it appears that the time t^=8.0215 days is a suspecious point of activity change. Fitting the data by two stage Omori-Utsu models, shows that the smaller AIC value, equal to -403.6618 is obtained for simple Omori-Utsu model. The results obtained are shown on Fig. 5. The results obtained for the two stage Omori-Utsu models are shown on Fig 5 (b), (c) and (d). The AIC values are include on the graphs.

Figure 5.

Graphs showing the cumulative number of aftershocks of Al Hoceima earthquake of 1994, with magnitude greater or equal to mc, fitted by a simple Modified Omori law – graph (a)- and two stage Modified Omori law, including secondary aftershock shown on the plot – graphs b, c and d – The parameters of each model and the AIC are given on the plot. (A) for AL Hoceima 1994 aftershock sequence and (B) for AL Hoceima 2004 afetsrhock sequence.

The analysis of the aftershock sequence of Al Hoceima 2004, give us that a suspecious point of aftershock activity change is given by t¯=2.402 days, the lower AIC equal to AIC= -3787.1548 is obtauined for the simple Omori-Utsu model as shown on the following Fig. 6(B). The results obtained for the two stage Omori-Utsu models are displayed on Fig. 5 (b), (c) and (d), including the AIC values. The principal reason of this, seen from the comparaison of the predicted and real cumulative curves would be that the secondary aftershock triggered by large aftershocks are not frequent enough in relation to their magnitude. It is well established that the residual analysis of point process is a "good" tool to evaluate the goodness of fit of the selected model to the data (Ogata, 1988; Ogata, 1999). The integral


of non-negative conditional intensity function produces a 1-1 transformation of the time from t to τ=Λ(t), so that the occurrence times t1,t2,t3,t4,t5,..................,tN are transformed into τ1,τ2,τ3,τ4,τ5,..................,τN. It is well known that τk,k=1,.....,N are distributed as a standard Poisson process. If λ(t|t) is the intensity function of the process actually generating the data, using the maximum likelihood conditional intensity λθ^(t|t), the corresponding τ^1,τ^2,τ^3,τ^4,τ^5,..................,τ^N called residual process (Ogata, 1988; 1999), provides a measure of teh deviation of the data from the hypothezed model. The residual analysis performed for the two aftershock sequences, trigged respectively by Al Hoceima earthquakes of 1994 and 2004 is shown on Fig. 6. The graphs display the Omori-Utsu model fit and the residual process for Al Hoceima 1994 series and for Al Hoceima 2004 series.

Figure 6.

Residual analysis to derive the best fitting for the two Al Hoceima aftershocks.

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.

Akaike Information Criteria
Model 1 Model 2 Model 3 Model 4 B. Model
Al Hoceima 1994 -639.1429 -609.7559 -609.276 -633.7404 Model 1
AlHoceima 2004 -3616.8231 -36020.7106 -3609.2007 -3610.6351 Model 1
El Asnam 1980 -47.1127 -32.118 -31.2841 -42.9912 Model 1
Zemouri 2003 -8973.1703 -9244.9968 -9342.0748 -9353.1883 Model 4
Laalam 2006 -293.5383 -274.8197 -275.9812 -288.698 Model 1

Table 2.

Akaike information criteria (AIC) of each model and for each aftershock sequence. Table gives in the last column the best model with the lower AIC.

Following figure 7.6 displays for El Asnam 1980, Zemouri 2003 and Laalam 2006 aftershock sequences the adjustment of the data with the appropriate model, as deduce from the Table2

Figure 7.

Adjustment of the cumulative number of events by the appropriate Omori-Utsu model, first stage Omori-Utsu model for (a) EL Asnam 1980 aftershock sequence and (b) for Laalam aftershock sequence. Two stage Omori-Utsu model for Zemouri 2003 aftershock sequence

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 t0 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


where ti is the time from the mainshock to the ith aftreshok, T the absolute temperature, R the universal gas constant, γ is a constant and S(ti) the cumulative stress drop. Since the seismic moment may be defined as M0=ΔσV, where V is the focal volume (Madariaga, 1979), the previous equation can be written as follow




Vk is the cumulative focal volume, V0 is the focal volume of the main shock, Vj is the focal volume of the jth aftershock and M0m and M0j are the seismic moment of the main shock and the jth 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 Ki 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:

(a) The validity of equation 22, which is adjusted to different aftershock data and plotted on Figure 8. The same figure shows the estimates a and b, their respective standard deviation σa and σb and the determination coefficient r2 of the linear adjustment on the semi-logarithmic scale. Solid curves are plots of the function y=a+bx where x=log(t). Dashed curves say (Δi), y=a+bx and (Δi+), y=a++b+x, correspond to the 95% confidence limits relative to the confidence intervals I(a)=[a,a+] and I(b)=[b,b+] of a and b, respectively. Namely, b+=b^+tn2(α/2)Sb^, b=b^tn2(α/2)Sb^; where tn2(α/2) is the α/2-quantile of the Student distribution with n2 degrees of freedom (α=0.05 in our case) with


ε^i=yiy^i are the regression errors and y^i=a^+b^xi the predicted yi values.

The confidence limite of a^ are calculated from those of b^ using the relation a^=y¯b^x¯.

(b) The validity of the definition of the constants a and b as expressed previously by the relation (23). The conditions (a) and (b) are analyzed first by the quality of the fit, shown on figure 8.

Figure 8.

Fit of aftershock sequence, the plot displays the 95% confidence limit of the regression, the parameters of the adjustment and the coefficient of the correlation are shown on the plot. For Al Hoceima series of 1994 and 2004 the results are shown on graph (a) and (b) respectively. The graph (c) displays the results of the El Asnam serie of 10 october 1980. Graphs (d) and (e) display the results for the series of Zemouri 2003 and Laalam 2006.

Taking into account the obtained coefficient of correlation r2 and given on each graph, the modele y=a+bLog(t) introduced by Marcellini (1997) fit very well our data.


5. 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 mms and its largest aftershock magnitude masmax is approximately constant, independently of the mainshock magnitude (Bath, 1965). That is,


Δmis typically around 1.2.

Extensive studies about the statistical variability of Δm have been carried out by several authors e.g Vere-Jones (1969), Kisskinger and Jones (1991) and Console et al. (2003). However, the law remains an open problem nowadays (Vere-Jones, 1969; Vere-Jones et al. 2005). In this study, we use the modified Bath's law proposed by Shcherbakov and Turcote (2004a). It is based on an extrapolation of the Gutenberg-Richter relationship for aftershocks. Namely, the magnitude of the largest aftershock consistent with the Gutenberg-Richter relationship for aftershocks is obtained by formally putting N(m*)=1 which yields to abm*=0. If Bath law is applicable to the inferres magnitude m*, the Gutenberg-Richter relationship takes the following form;




combining this last relation with the empirical relation of the energy dissipated


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)




This relation is applied directly to describe the link between the energy radiated by the mainshock Ems and the moment magnitude of the mainshock mms,


following Shcherbakov et al.(2004a), the total energy radiated during the aftershock sequence Eas is obtained by integrating over the distribution of aftershock till the inferred largest magnitude as upper bound in the integration. The upper bound of the integration is the largest aftershock magnitude unferred from the Gutenberg-Richter relationship and is denoted m*. The total radiated energy in the aftershock sequence is obtained by integrating over the distribution of aftershock as,


m* is the largest aftershock magnitude inferred from the Gutenberg-Richter law. Taking the derivative of the modified Bath law,


by combining the former equations, Eq. 33 and 34, we obtain


in addition, different version of the above equation, Eq 35 is obtained if we use the equation giving E(m) as a function of E0 and b value, we get


taking into account the modified Bath's law, we obtain


the ratio of the total radiated energy by the aftershocks Eas to the total energy radiated by the mainshock Ems 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.

Ratio of the elastic energy released
b±σb Δm* EasEas+Eas
Al Hoceima 1994 1.07 ± 0.07 1.10 0.05
Al Hoceima 2004 1.13 ± 0.05 0.50 0.35
El Asnam 1980 0.82 ±0.10 1.20 0.02
Zemouri 2003 1.10 ±0.04 0.82 0.14
Laalam 2006 0.99 ±0.10 1.70 0.01

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.


6. 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 D2. This D2 - value is an extension of the Euclidian dimension and measures the degree of clustering of earthquakes. In the two dimensional space, D2 can be a decimal number and ranges from 0 to 2.0. Therefore, the distributions characterized by different fractal dimensions defiine different clustering of events in space. In the two dimensional space, as D2 approaches 1, the distribution of events approaches a line (Euclidean dimension equal to 1). The same occurs for the distribution of events along a fault. Alternatively, as D2 approaches 2, the distribution of events tends to be uniform on the plane (Euclidean dimension equal to 2). When the D2-values approaches 0, the distribution is concentrated in a single point (Beauval et al. 2006; Spada et al. 2011). In this study, the fractal dimension is estimated using the correlation dimension (Grassberger and Procaccia, 1983)


where r is the radius of the sphere considered in the study, and C(r) is the correlation integral given by


where N is the number of points in the analysis window, xi; i=1,...., N are the coordinates of the epicenters, and H(.) the Heaviside step function, H(x)=0forx0,

H(x)=1forx0. The quantity C(r) is equivalent to the probability that two points will be separated by a distance less than r. The correlation integral is theoretically proportional to the power of D2, i.e C(r)rD2, where D2 is the correlation dimension equals the second generalized Renyi dilmension d2 (Molchan and Kronod, 2009). To estimate D2 from the correlation integral (Spada et al. 2011), we plot C(r) versus r on the log-log scale, Fig. 6., then we use the least square method to fit the data over the region where C(r) is linear, which corresponds to the gradient of straight line of the resulting plot of log10C(r) against log10(r).

Figure 9.

Graph displaying the plot of the correlation integrale in log-log scale, with the 95% of confidence limit (in dashed lines), the graphs show also, the plot the slope, related to the first derivative of the correlation integral for Al Hoceima 1994 and 2004 aftershock sequences.

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, rd and rs 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 rd/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

Figure 10.

Graph showing the correlation integral Log10C(r)vsLog10(r).

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 Sp is the slip over primary fault and S represents the total slip over the fault-system.

The obtained results are shown on Table 3, we deduce that during Al Hoceima earthquake of 1994, 62% of the total slip accomodates the primary rupture, 60 % during Al Hoceima 2004 earthquake. During the El Asnam earthquake of 10 October 1980, 59% of the total slip accomodates the primary fault segment, during Zemouri earthquake of 2003 this ratio has been estimate to 57% and 45% during Laalam earthqiuake of 2006. It is important to point out that in each case the remainder of the slip is distributed over the secondary rupture.

Fractale dimension and ratio of the slip
D2±σD2 Range Sp/S
Al Hoceima 1994 1.60 ± 0.05 1.79 - 12.56 0.62
Al Hoceima 2004 1.67 ± 0.04 2.16 - 17.27 0.60
El Asnam 1980 1.70 ±0.09 12.90 - 31.35 0.59
Zemouri 2003 1.79 ±0.02 1.60 - 10.00 0.57
Laalam 2006 2.13 ±0.02 1.57 - 8.90 0.45

Table 4.

Fractale dimension D2 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 x1,x2,.........,xn 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 x1,x2,.........,xn; the quadratic uncertainty is then given by


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




assuming that the minimum of M0(h) is in the vicinity of the minimum of E[M0(h)] and then it is in the vicinity of E{+[f^(x)f(x)]2dx}. Also, for n large, (n1)n; finaly the problem consites to minimize the simplified score function defined as follow, we replace M0(h) by M1(h) given by


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 f^i(.,h) is the kernel density calculated from the sample x1,x2,...........,xi1,xi+1,........,xn. We use the Silverman multi-modal tests to estimate the number of true bumps. This test is as follow T(k) test the null hypothesis H0k: " f has k bumps " against the alternative H1k: " f has more than k bumps". Under the hypothesis H0k, the smoothing parameter of the Gaussian kernel density estimate f^(.,h) is given by


The test T(k) is inspired from the fact that big values of parameter hcrit(k) reject the null hypothesis H0k. The following therorem by Silverman (1981) gives a characterisation

Theorem (Silverman, 1981)

The kernel density estimate f^(.,h) with Gaussian kernel, has more than k bumps if and only h<hcrit(k).

The test T(k) is constructed by simulating N statistics hcrit(1)(k);hcrit(2)(k);hcrit(3)(k);.........;hcrit(N)(k) from N smoothed bootstrap samples of size n obtained from x1,x2,.........,xn.

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 Efron formula,


where (xI(i);i=1,n) is a bootstrap sample of size n simulated from the sample x_=(x1,x2,.........,xn); x¯ and σ2 are the mean value and variance of the sample x_, and (εi;i=1,n) is a randomly simulated sample from the standard normal distribution.

if hcrit*(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 T(1),T(2),T(3),........., in this order until the null hypothesis is not rejected. We used the previous methods, to derive the distribution of the inter-event distances. The results obtained for the aftershock sequence triggered by the Al Hoceima 1994 earthquake are shown on Fig. 11.

Figure 11.

Kernel density estimated for Al Hoceima 1994 aftershock sequence (a) using the rule of thumb (b)using cross validation smoothing and (c) Estimated density obtained with one mode critical bandwidth hcrit=6.4

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(b). The estimated density obtained with one mode critical bandwidth equal to hcrit=6.4 is shown on Fig. 11(c). On the three Figures a concentration of the inter-event distances around the value 8.9 km, which correspond to the mode of the estimated distribution of probability. The results obtained for the Al Hoceima 2004 aftershock sequence are shown on the Fig. 12. The estimated density using the rule of thumb shown on Fig 12(a) displays a concentration of the inter-event distances around the value 9.45 km, the estimated density using the cross validation optimal smoothing has been obtained using a bandwidth parameter h=1.0. The density obtained displays a concentration of the inetr-event distances around the value 7.5 km. Although, the estimated density with one mode critical bandwidth obtained equal to hcrit=11, shows a concentration of the inter-event distances around the value 14.14 km.

Figure 12.

Kernel density estimated for Al Hoceima 2004 aftershock sequence (a) using the rule of thumb (b)using cross validation smoothing and (c) Estimated density obtained with one mode critical bandwidth hcrit=11

For the aftershock sequence triggered by the 21 May 2003 Zemouri earthquake (Mw 6.9), the obtained results are shown on Fig. 13.

Figure 13.

Kernel density estimated for Zemouri 2003 aftershock sequence (a) using the rule of thumb. (b) using cross validation smoothing and (c) Estimated density obtained with one mode critical bandwidth hcrit=6.6

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 hcrit=6.6; displays a concentration of the inter-event distances around 9.0 km.

Figure 14.

Kernel density estimate for Laalam 2006 aftershock sequence ((a), (b) and (c)). The Graphs on (d) and (e) display the kernel density estimate for El Asnam 1980 aftershock sequence.


7. 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 mc 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 mc, we used the maximum likelihood to estimate above the threshold completeness mc. 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.8mbLg. 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 D2 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.


  1. 1. AkaikeH1974A new look at the statistical model identification, IEEE Trans. Autom. Control 19716723
  2. 2. AkiK1965Maximum likelihood estimate of b in teh formula Log N = a- bM and its confidence limits,. Bull. Earthq. Res. Inst. Tokyo Univ. 43237239
  3. 3. ArmoreseD2007Applying a change of point detection method on frequency-magnitude distributions.Bull. Seismology Society of America, doi.
  4. 4. BathM1965Lateral inhomogeneties in the upper mantle. Tectonophysics, 2483514
  5. 5. BeauvalCSHainzland FScherbaum2006The Impact of the spatial uniform distribution of seismicity on probabilistic seismic-hazard estimation. Bulletin Seismology Society of America. 96n6 24652471
  6. 6. BeldjoudiHAGuemacheAKherroubiFSemmaneK. AYelles-chaoucheHDjellitAAmraniand AHaned2009The Laalam (Bejaia, Nort-East Algeria) Moderate earthquake (Mw=5.2) on March 20, 2006. Pure and Applied Geoph. DOIs00024-009-0462-9.
  7. 7. BenderB1983Maximum-Likelihood estimation of b-values for magnitude grouped data, Bull. Seismol. Soc. Am. 73, n°3831851
  8. 8. CalvertAGomezFSeberDBaranzagiMJabourNIbenbrahimAand ADemnati1997An integrated geophysical investigation of recent seismicity in the Al Hoceima region of North Morocco. Bulletin Seism. Soc. of Am. 87. 637651
  9. 9. ConsoleRA. MLombardiMMurruand RRhodes2003Bath’s law and the self-similarity. J. Geophys. Res. 108, 2128.
  10. 10. Daley and Vere-Jones2003An introduction to the theory of point process, 2nd Ed., 1Springer-verlag. New-York. Berlin Heidelberg.
  11. 11. El AlamiS. OTadiliB. ACherkaouiT. EMedinaFRamdaniMAit-brahimLand MHanafi1998The Al Hoceima earthquake of May 26, 1994 and its aftershocks : a seismotectonic study. Annali di Geofisica, 41519537
  12. 12. EnevaM1996Effect of limited data sets in evaluating the scaling properties of spatially distributed data : an example from minning-induced seismic activity, Geophys., J. Int., 124773786
  13. 13. EnescuBand KIto2002Spatial analysis of the frequency-magnitude distribution and decay rate of aftershock activity of the 2000 Western Tottori earthquake. Earth Planets Space 54847859
  14. 14. EnescuBJMoriMMasatoshiand YKano2009Omori-Utsu law c-values associated with recent moderate earthquakes in Japan, Bull. Seismol. Soc. of Am. 99, 2A, 884891
  15. 15. FrohlichCand S. DDavis1990Single-link cluster analysis as a method to evaluate spatial and temporalproperties of earthquake catalogues. Geophy. J. Inter. 1001932
  16. 16. Galindo-zaldivarJAChalouanOAzzouzCSanDe GaldeanoFAnahnahLAmezaPRuanoAPedreraARuiz-constanCMarin-lechadoMBenmakhloufA. CLopez-garridoMAhmamouRSajiF. JRoldan-garciaMAkli., and A. Chabli. 2009Are seismological and geological observations of the Al Hoceima (Morocco Rif) 2004 earthquake (M=6.3) contradictory. Tectonophysics,
  17. 17. GardnerJ. Kand LKnopoff1974Is the sequence of earthquake in southern California, with aftershock removed, Poissonian ? Bull. Seism. Soc. Am., 64 (5), 1363-1367.
  18. 18. GoltzC1998Fractal and chaotic properties of earthquake, in Lecture Notes in Earth Sciences, Springer, New York, 175 pp.
  19. 19. GrassbergerPand IProcaccia1983Measuring the strangeness of strange attractors. Physics D, 9189208
  20. 20. GutenbergRand C. FRichter1944Frequency of Earthquake in California. Bull. Seismol. Soc. Am. 34. 158188
  21. 21. GuoZand YOgata1997Statistical relation between the parameters of aftershocks in time, space and magnitude. Journal of Geophys. Res. 102(B2).28572873
  22. 22. GuttorpPand DHopkins1986On estimating varying b-values, Bull. Seismol. Soc. Am. 76 n°3889895
  23. 23. HamdacheMPelaézJ. Aand K. Yelles Chaouche. 2004The Algiers, Algeria earthquake (Mw 6.8) of the 21 May 2003: Preliminary report. Seism. Res. Letters, 75n3. May/June 2004.
  24. 24. HamdacheMPelaézJ. ATalbiAand López Casado, C. 2010A unified catalog of main earthquakes for Northern Algeria from A.D. 856 to 2008. Seism. Res. Lett. 81732739
  25. 25. KaganY. Y2004Short-term proprieties of earthquake catalogs and models of earthquake source. Bull. Seismol. Soc. Am. 94 412071228
  26. 26. KhattriK. N1995Fractal description of seismicity of India and inferences regarding earthquake hazard. Curr. Sci., 69. 361366
  27. 27. KisslingerCand L. MJones1991Proprieties of Aftershocks in Southern California. J. Geophy. Res. 1032445324
  28. 28. MaoucheSMMeghraouiCMorhangeSBelabbesYBouhadadand HHaddoum2011Active coastal thrusting and folding, and uplift rate of the Sahel Anticline and Zemouri earthquake area (Tell Atlas, Algeria). Tectonophysics, 509, 6980
  29. 29. MarcelliniA1995Arrhenius behavior of aftershock sequences. J. Geophys. Res. 10064636468
  30. 30. MarcelliniA1997Physical model of aftershock temporal behavior. Tectonophysics 277137146
  31. 31. MarzocchiWand LSandri2003AReviewand new insights on the estimation of the b-value and its uncertainty. Annals of Geophysics. 46n6.
  32. 32. MikumoTand TMiyatake1979Earthquake sequences on a frictional fault model with non-uniform strenghs and relaxation times, Geophys. Journal R. Astron. Soc., 59497522
  33. 33. MogiK1962Study of elastic shocks caused by the fracture of heterogeneous materials and its relation to the earthquake phenomena, Bull., Earthquake Res., Inst., Univ., Tokyo, 40125173
  34. 34. MolchanGand TKronod2009The fractal description of seismicity, Geophs. J. Int. 179. n 317871799doij.1365.
  35. 35. NerenbergM. A. Hand CEssex1990Correlation dimension and systematic geometric effects, Phys. Rev. A. 42, 70657074
  36. 36. NocquetJ. Mand ECalais2003Crustal velocity field of Western Europe from permanent GPS array solutions, 1996-2001. Geoph. J. International, 1547288
  37. 37. NyffenggerPand CFrolich1998Recommandations for determining p values for aftershock sequence and catalogs. Bull. Seismol, Soc. Am. 88n0 5 11441154
  38. 38. NyffenggerPand CFrolich2000Aftershock occurrence rate decay properties for intermediate and deep earthquake sequences. Geoph. Res. Lett. 2712151218
  39. 39. OgataY1983Estimation of the parameters in the modified Omori formula for aftershock frequencies by the maximum likelihood procedure. J. Phys. Earth. 31115124
  40. 40. OgataY1992Detection of precursory relative quiescence before great earthquake through a statistical model. Geophys. Res. 971984519
  41. 41. OgataYand KKatsura1993Analysis of temporal and spatial heterogeinity of magnitude frequencey distribution inferred from earthquake catalogue. Geophys. J. Int. 113, 727738
  42. 42. OgataY1999Seismicity Analysis through Point-Process Modelling : A Review. Pure and Applied Geophysics, 155. 471507
  43. 43. OgataYL. MJonesand SToda2003When and where the aftershock activity was depressed : Contrasting decay patterns of the `proximate large earthquake in southern California. J. of Geophysics Research, 108n0 B6 2318, doiJB002009.
  44. 44. OlssenR1999An estimation of the maximum b values in the Gutenberg-Richter relation. Geodynamics, 27547552
  45. 45. OmoriF1894On the aftershocks of earthquake. J. Coll. Sci. Imp. Univ. Tokyo, 7. 111120
  46. 46. OuyedMMMeghraouiACisternasADeschampADorelFFrechetRGaulonDHatzfeldand HPhilip1981Nature, 29258182631
  47. 47. PeláezJ. AChourakMTadiliB. ABrahimL. AHamdacheMLópez Casado, C., and Martínez Solares, J.M. 2007A catalog of main Moroccan earthquakes from 1045 to 2005. Seismological Research Letters 78614621
  48. 48. Peláez, J.A., M. Hamdache and Sanz De Galdeano. 2012. A spatially smoothed seismicity forecasting model for Mw5.0earthquakes in northern Algeria and Morocco. 15 World Conf. Earth. Eng. 24 - 28 September 2012, Lisboa, Portugal
  49. 49. ReasenbergP1985Second-order moment of central California seismicity, 1969-82, J. Geophys. Res., 9054795495
  50. 50. RydelekP. Aand I. SSacks1989Testing the completeness of earthquake catalogues and the hypothesis of self-similarity, Nature, 337251253
  51. 51. SandriLand WMarzocchi2005A technical note on the bias in the estimation of the b-value and its uncertainty through the least squares technique. Annals of Geophysics.
  52. 52. Sanz de GaldeanoC., 1990Geologic evolution of the Betic cordilleras in the Western Mediterranean, Miocene to Present. Tectonophysics 172107119
  53. 53. ShcherbakovRand D. LTurcotte2004aA modified form of Bath’s law. Bull. Seismol. Soc. Am. 9419681975
  54. 54. ShcherbakovRD. LTurcotteand J. ERundle2004bA generalized Omori´s law for earthquake aftershock decay. Geophy. Res. Lett. 31. L11613, doiGL019808.
  55. 55. ShcherbakovRD. LTurcotteand J. ERundle2005Aftershocks Statistics. Pure and Applied Geophysics 16210511076
  56. 56. ShcherbakovRand D. LTurcotte2006Scaling properties of the Park-field aftershock sequence. Bull. Seismol, Soc. Am. 94. S376S384.
  57. 57. ShiYand BBolt1982The standard error of the magnitude-frequency b value. Bull. Seism. Soc. Am. 72n°5 16771687
  58. 58. SilvermanB. W1986Density estimation for statistics and data analysis. Chapman and Hall, London.
  59. 59. SpadaMSWeimerand EKissling2010Quantifying a potential bias in probabilistic seismic hazard assessment: Seismotectonic zonation with fractal properties. Bull. Seism. Soc. Am. 101n 6 26942711
  60. 60. StichDFMancillaDBaumontand JMorales2005Source analysis of the Mw 6.3 2004 Al Hoceima earthquake (Morocco) using regional apparent source time functions. J. Geoph. Res. 110 (B06306) doiJB003366.
  61. 61. TahirMJ. RGrassoand DAmorèse2012The largest aftershock: How strong, how far away, how delayed?. Geophysical Research letters, 39L04301, doi:GL050604, 2012.
  62. 62. ThatcherWand T. CHanks1973Source parameters of southern California earthquake. J. Geophys. Res. 7885478576
  63. 63. TintiSand FMulargia1987Confidence intervals of b-values for gropuped magnitudes. Bull. Seismol. Soc. Am. 77, n°621252134
  64. 64. TurcotteD1997Fractals and chaos in Geology and Geophysics, Cambridge University Press, New York, 416 pp
  65. 65. UtsuT1961A statistical study on the occurrence of aftershocks. Geophysics, 30. 521605
  66. 66. UtsuT1965A method for determining the value of b in a formula log n = a- bM showing the magnitude-frequency relation for earthquakles, Geophys, Bull. Hokkaido Univ. 1399103
  67. 67. UtsuT1969Aftershocks and earhquake statistics (I)- Some Parameters which characterize an Aftershock Sequence and their Interaction. Journal Fac., Sci., Hokaido, Univ., Ser., VII (Geophys.), 3129195
  68. 68. UtsuT1971Aftershocks and Earthquake Statistics (III). Analyses of teh distribution of earthquakes in magnitude, Time and Space with special consideration to clustering characteristics of earthquake occurrence. Jpurnal of the Faculty of Science, Hokkaido Univ. Ser. VII, Geophysics, Vol. III, n°5.
  69. 69. UtsuTYOgataand R. SMatsura. 1995The centenary of the Omori formula for a decay law of aftershock activity. J. Phys. Earth 43. 133
  70. 70. Vere-jonesD1969A note on the statistical interpretation of Bath`s law. Bulletin Seismol. Soc. of Am. 69n4, 15351541
  71. 71. Vere-jonesDJMurakamiand AChristophersen2005A further note on Bath’s law, The 4th International Workshop on Statistical Seismology.Tokyo. Japan.
  72. 72. VidalF1986Sismotectónica de la región Béticas-Mar de Alborán. PhD Thesis. Universidad de Granada. 457 p.
  73. 73. WiemerSand R. FZuniga1994Zmap. A Software package to analyse seismicity (abstract), EOS Trans. AGU 75 (43), Fall Meet. Suppl., 456.
  74. 74. WiemerSand KKatsumata1999Spatial variability of seismicity parameters in aftershock zones. J. Geophys. Res. 104 (B6), 13, 135-151.
  75. 75. WiemerSand MWyss2000Minimum magnitude of completness in earthquake catalogs : Examples from Alaska, the Western United States, and Japan. Bull. Seismol. Soc. Am. 90859869
  76. 76. WiemerS2001A software package to analyze seismicity: ZMAP, Seismol. Res. Lett., 72374383
  77. 77. WiemerSand MBaer2000Mapping and removing quarry blast events from seismicity catalogs, Bull. Seism. Soc. Am., 90525530
  78. 78. WiemerSand MWyss2000Minimum magnitude of complete reporting in earthquake catalogs: examples from Alaska, the Western United States, and Japan, Bull. Seism. Soc. Am. 90859869
  79. 79. WoessnerJand SWiemer2005Assessing the Quality of Earthquake Catalogues: Estimating the Magnitude of Completeness and Its Uncertainty, Bull. Seismol. Soc. Am., 95684698doi:10.1785/0120040007,

Written By

M. Hamdache, J.A. Pelàez and A. Talbi

Submitted: 04 May 2012 Published: 20 March 2013