Conventional and Fractal Variogram Based on Time—Space Analysis of Seismicity Distribution—Case Study: Algeria Seismicity (1673-2010)

Moh-Amokrane Aitouche^*, Mabrouk Djeddi

Department of Geophysics, Faculty of Hydrocarbons and Chemistry, Earth Physics Laboratory, University of Boumerdes, Boumerdes, Algeria.

DOI: 10.4236/gep.2018.611012                

Abstract

Geostatistics belongs to the wide class of statistical methods. It is used and applied to analyze and predict the values associated with spatial or spatio-temporal phenomena such as seismological events. Thus, in addition to its adaptability to perform spatial data analysis based on the principles of variography, many other methods more representative of the spatial distributions of the epicentres have been experimented. The unidirectional or isotropic fractal variogram associates the criterion of variance with the fractal dimension which is a better descriptive parameter of the spatial organization of events. The same analysis procedure carried out through a directional or azimuthal variogram introduces the context of a preferential direction of earthquake evolution. Moreover, b-value is a privileged seismological parameter that can be studied alone or in combination with other more advanced geostatistical analysis factors such as fractals, fractal dimensions and the anisotropic variogram. All these concepts were used as methodology for a protocol of analysis of the catalog of the Algeria seismicity including 1919 events for the period extending from 1673 to 2010. This catalog is provided by the compilation of partial catalogs synthesized by the CRAAG (Research Centre in Astronomy, Astrophysics and Geophysics).

Keywords

Algeria, Seismicity, Catalog, b-Value, Geostatistics, Semiv Variogr, Am Fractal

1. Introduction

Many natural phenomena are identifiable by measurements or recordings. Thus, the seismological data can be processed using stochastic approach or on the basis of existing uncertainties. More generally, statistical methods like Geostatistics are proving to be highly adapted tools to make it possible to provide an efficient assessment of the analysis parameters. The seismological data recorded in a catalog are identifiable through three types of dimensions: the space dimension (i.e. longitude and latitude) the magnitude and the time of occurrence (Wierner, 1996; Ogata & Katsura, 1993) . In the strictly geostatistical context, these data are assimilated to time series obeying the laws of distribution, clustering and scale invariance. The seismological data recorded in a catalog are identifiable through three types of dimensions, a dimension of space (longitude and latitude), magnitude and time of occurrence (Wiemer & Wiss, 2002) . On the purely geostatistical level, these data are assimilated to time series whose distribution laws, clustering mode and scale invariance status will be determined in the fractal case. These characteristics are elucidated only in relation to the parameters of analysis imposed by the objectives of the study. Various statistical techniques for studying local or regional seismicity have been applied to catalog data across a restricted geographic area or country. The frequency-magnitude distribution provided by the associated Gutenberg-Richter relationship (Gutenberg & Richter, 1942) and the fractal dimension concept (Hirata, 1989; Legrand, 2002) is considered as a convenient and objective approach in the seismicity analysis.

Conventional statistical analysis provides only a limited set of tools because it could be unaware certain no obvious aspects of the data to be analyzed (Cressie, 1993) . So the French mathematician Georges Matheron introduced geostatistics as a set of probabilistic methods designed to understand the structure of variables regionalized and distributed in a real space like seismicity data (Mathéron, 1970; Schaefer et al., 2014) .

Contrary to conventional statistics, methods, geostatistics takes care through adapted mechanisms of study and specific concepts describing the complexity and irregularity of the phenomena (Goovaerts, 1997; Mathéron, 1970) . Thus, the more recent statistical approaches are variogram (Olea, 2009) and directional or azimuthal variogram. Fractal structures, fractal variogram and the principle of spatial heterogeneity are associated with the classical concepts of variance, correlation, and random distribution. All these arguments require the implementation of specific methods combining geostatistical processing, variography and fractal analysis. The materials and methods included in this study aim to demystify the seismic catalog by highlighting the intrinsic parameters of earthquakes. The number of seismic events extracted from the catalog before declustering is around 1919 events reduced to 1879 after eliminating the non-significant smallest events and after the declustering step. Because of their direct connection with sismotectonic activity, the b-value extracted from the Gutenberg-Richter relationship, the variogram representation, the fractal dimensions and the interconnections between these parameters will constitute the framework of this study.

2. Regional Geodynamical Context

Mediterranean Sea belongs to an active and permanent geodynamic context resulting from the collision of the two African and Eurasian plates inducing a frequent but moderate seismicity. Indeed, the basins delimiting the extent of the Mediterranean form a perimeter zone of separation of these two geodynamic entities. This limit is characterized by a strong geological complexity induced by a sustainable and diversified geodynamics. This geodynamics includes alternating sequences of extensions and compressions. The extension phases lead to the appearance of new oceans whereas the compression phases led to the appearance of subduction zones and to the creation of collision chains. This geological complexity continues until recently in terms of deformation and seismotectonic activity. Overall, these phenomena are local or regional but exert an impact on the deformation field related to the convergence of African and Eurasian tectonic plates (Figure 1).

The coexistence of the compressive and distensile deformation regimes is the driving force of the NW direction convergence of the two plates. However, the speed of convergence is not uniform; it is estimated at between 3 and 8 mm/yr in the eastern Mediterranean and at 2 and 5 mm/yr in the western Mediterranean (Figure 2).

3. Seismotectonic Framework of Algeria

Algeria country is a part of the north-west of the African continent. It also belongs to the eastern Ibero-Maghrebian region. In the regional tectonic context, Algeria is located on the edge of the African and Eurasian plates, and it undergoes the effects induced by the convergence dynamics of these two plates (Kenzi, 1972) . The seismicity of Algeria is located on the northern fringe of the country formed by four morpho-structural domains that are the Tellian Atlas, the High Plateaus, the Saharan Atlas and the northern part of the Saharan Atlas platform (Figure 3). Only the first two domains cited feel the effects of the inter-plate movements confirming their high seismicity whereas the Algerian south is relatively spared by the seismic activity.

Figure 1. Convergence of African and Eurasian tectonic plates-zone of seismicity.

Figure 2. Geodynamical setting of western Mediterranean Sea (after Carminati et al. 2012) .

Figure 3. Main Geographic regions of Algeria.

The origin of this proven seismicity comes from the compressive regime to which the Tell is subjected, with a shortening of direction NNW-SSE. Earthquakes are generally generated by folds, faulty folds and especially by NE-SW oriented faults. The frequency and magnitude of the seismicity are important on the Tellian Atlas while this seismicity is ostensibly attenuated towards the southern Algeria. Moderate earthquakes are recorded in the Saharan Atlas and in the part of the Saharan platform. The frequency and magnitude of the seismicity are more important on the Tellian Atlas but this seismicity is disparate towards the southern Algeria region. However, historically, some moderate earthquakes have been recorded the Saharan platform.

4. Identification of Earthquakes in Algeria

Northern Algeria is characterized throughout its extent by seismicity caused by the confrontation of the African and Eurasian Plates. This seismicity involves the occurrence of weak to moderate earthquakes although sometimes violent earthquakes can occur. These earthquakes generate many disasters like the earthquakes of El-Asnam of 10/10/1980 (Benouar, 1994) and more recently the May 21, 2003 Boumerdes earthquake (Aitouche et al., 2012) . On a larger time scale, a large number of earthquakes were recorded, including those used in the catalog on which this study was based for a period ranging from 1673 to 2010 (Figure 4). One can cite for this purpose, the earthquakes of Setif (419), supposed to be the first known historical earthquake, then Algiers (1365 and 1716), Oran in 1790, and Gouraya in 1891. In a more recent period, we can cite, for example, the Orlénvilles earthquakes of 9/9/1954, El-Asnam (actually Ech-Cheliff) of 10/10/1980, Constantine of 27/10/1985.

Figure 4. Spatial distribution of earthquakes―Primary statistical settings.

5. Seismic Catalog Compilation and Selective Sorting of Data

It is not easy to produce a homogeneous and reliable seismic catalog over such a long period of time (1673-2010). The catalog used in this study is available from the Centre for Research in Astronomy, Astrophysics and Geophysics (CRAAG―Algiers). It is the result of a synthesis and a compilation of partial catalogs each attached to a precise interval of time (Benouar, 1994) . In addition to the catalog provided by the CRAAG, several authors contributed to the final version of the seismic catalog as used in this study. So, a number of 1919 seismic events before the declustering sequence and covering the period 1673-2010 for the northern region of Algeria (longitude (Longitude 2.85˚W - 10˚E; latitude 31.28˚S - 37.80˚N) defines the core of the earthquake catalog. The main obstacle to the construction of the catalog lies in an objective conversion of intensities into surface magnitudes.

Figure 5 is a perspective representation of the distribution of earthquakes as a function of longitude, latitude and magnitude. The panoramic view of the events illustrates perfectly the concentration of the seismic activity along the Algerian Tell with a dominant magnitude of 7.5 as shown by the major peak of this distribution.

6. Seismicity and Tectonic Setting of Northern Algeria

Northern Algeria is affected all along its Tellian part by a recurring seismicity but diversified in terms of periods of occurrence and magnitude (Bellalem et al., 2008) . Local or regional tectonics materialized by a network of active faults is the main source of this recurrent seismicity. Figure 6 displays the distribution of epicentres of the seismic events and the corresponding seismal curves. Some of these epicentres present alignments of orientations NW-SE or NE-SW identifiable with the tectonic regime of northern Algeria. The isoseismal curves show well localized accumulation points, thus reflecting a strong local seismicity. The

Figure 5. 3D-Perspective view of the spatial earthquake distribution.

Figure 6. Isoseismal curves―spatial positioning of epicentres.

circled point has an anomaly that differs from all other points. It represents an area with high seismic potential. It is also worth noting the isolation of some epicentres which could indicate a sporadic seismicity due, perhaps, to a fault reactivation. Areas with low seismicity are either not visited by the isoseismal curves or by their relatively wide spacing.

7. Catalog Declustering

The occurrence of any global seismic event is generally composed of two specific partial events including the major seismic event followed by a sequence of induced earthquakes also known as aftershocks. Under the assumption that the seismological processes generating these two phenomena are statistically Gaussian, they are then supposed to be independent of each other with a random spatial distribute. It then becomes necessary to remove these aftershocks classified as noise. The declustering of the catalog is the first step to perform before any processing operation. The declustering is conceived as a process of separating the catalog data into foreshocks, mainshocks and aftershock events. There are several algorithms for removing the aftershocks like Rosenberg algorithm (1975). Currently, the algorithm proposed by Gardner and Knopoff (1974) is widely used due to its availability and simplicity. This algorithm is based on an empirical identification of the aftershock sequence using time and space windows that vary as a function of the mainshock magnitude.

7.1. Windowing Methods

The Gardner-Knopoff algorithm is based on the windowing approach to identify the mainshock and the aftershocks. For each earthquake listed in the catalog with a magnitude M the associated seismic events are identified as aftershocks if they belong to a specified time interval and within a range of corresponding distance. Thereafter, the time windows will be reset according to the magnitude value of the largest magnitude in the sequence. An approximation of the size of the temporal window T and distance D according to the Gardner and Knopoff algorithm is given by the following equations

$T\left(\text{Days}\right)=\{\begin{array}{l}{10}^{0.032M+2.7389}\text{if}M\ge 6.5\\ {10}^{0.5409M-0.547}\text{else}\end{array}$ (1)

$D\left(\text{Kms}\right)={10}^{0.1238M+0.983}$ (2)

7.2. Declustering of the Algerian Seismic Catalog

The Algerian seismicity catalog was declustered using Zmap software (Wierner, 2001) running under Matalb language (Wierner, 1996) . The Gardner-Knopoff’s algorithm is implemented with the possibility to vary the declustering parameters. This operation provided the map shown in Figure 7 which displays the seismic events of the original catalog after declustering as well as the parameters used to draw this map. We notice that only 20 clusters were built by the software containing in all 40 seismic events deduced from the 1919 initial number of earthquakes. The study will therefore be carried out using the difference i.e. 1879 earthquakes. This study was also carried out by the use of other software for the geostatistical data processing. Among them, the Zmap package running under Matlab language [28] and the GS+ (Geostatistics for the Environment) software.

8. First Results of the Geostatistical Analysis

The north of Algeria is characterized practically throughout the Tellian part by intense seismicity induced by the collision of the African and Eurasian plates. This seismicity generally results in moderate to weak earthquakes, although severe strong earthquakes may occur. An abnormal high magnitude ( $M_w=7.5$ ) was cited in the catalog. Figure 2 and Figure 3 display the spatial distribution of

Figure 7. Algeria catalog of seismicity after declustering.

epicentres of all earthquakes that occurred from the year 1673 to the year 2010. They illustrate a scattering of the epicentres along the entire Tellian band.

The cumulative number vs. time curve as shown in the figure (small squares) depicts a regular increase in the number of earthquakes for the period considered (1673-2010) with two remarks, however. The first relates to the presence of irregularities comparable to fluctuations in the number of earthquakes in the interval 1850-1950 (figure (small squares) probably due to a regional relaxation of the stress field. This state may be correlated in the same time with the frequency histogram (Figure 8) and the time histograms (Figure 4(c)). The earthquakes in Algeria are not uniformly distributed over time. At seismically quiet periods, phases of strong seismic activity can be followed. This means that a great seismic activity would have occurred. The earthquakes in Algeria are shallow. Their deep is estimated at a few tens of kilometres which make them devastating in case of occurrence (Figure 9).

For example, the destructive earthquake of magnitude $M_w=6.8$ which hit the balneal town earthquake of Boumerdes (east of the capital Algiers) on May 21, 2003 is ranked among the strongest seismic events that occur in this region. It was a shallow quake (10 km in depth) which caused important material and human damage in the epicentral region. Shallow earthquakes are caused by fault breaks or continental inter-plate movements. This may explain the large spread and the greater damage caused by some quakes that occurred in the Algerian Tellian region.

9. Frequency-Magnitude Distribution

The relationship between magnitude and frequency of occurrence of past earthquakes has been described by the Gutenberg-Richter power law (Gutenberg & Richter, 1942) as given by n the following equation

Figure 8. Magnitude variation vs. time.

Figure 9. Shallow distributions of quakes.

${\log }_{10}N\left(>M\right)=a-bM$ (3)

where $N\left(>M\right)$ indicates the cumulative number of earthquakes that have a magnitude greater than M while a and b are specific real positive coefficients which vary in a fixed time and space window (Crownover, 1995) . In a seismotectonic context, the constant a is related to the level of seismicity and to the seismic activity while the parameter b also known as b-value describes the number of the large seismic events that occurred or indirectly the level of the accumulated stress. b-value may be also used for interpretation of geological structures and to understand the faulting systems. The changes of b-value are an indicator of the level of stress in the area of study. Globally the b-value is close to 1.0, but experimentally it can vary from 0.5 to 1.5. Deviations of this parameter are then frequently observed and a mapping of the b-value variations can help with a better seismological investigation. Graphically, the b-value is calculated from the slope of the line representing the Equation (3) as showed in Figure 10. Applied to the catalog of Algeria’ seismicity, the slope of the red line corresponds to $b=0.773$ +/− 0.06 i.e. $0.713\le b\le 0.833$ . These values of b thus deviate from the standard value 1.0 with a tendency towards the low values. This is an indication of stress accumulation. The earthquakes in Algeria are shallow. Their deep is estimated at a few tens of kilometres which make them devastating in case of occurrence (Figure 9).

Figure 10. Cumulative numbers vs. magnitude.

10. Spatio-Temporal Variations of the b-Value across the Algerian Tell

Since b-value is related to the size of the seismic event distribution, it can thus define a measure of the frequency of weak or strong earthquakes that occurred in the region under study (Wackermagel, 2003; Mc Nutt et al., 2004) . Moreover, b-value is also an indicator of the tectonic activity and consequently of the field of stress. Figure 11 shows how the b-value varies with time for the considered catalog. The grey area highlights the zone of decrease of the b-value which corresponds to the time interval 1850-1950. It is precisely the same interval observed in Figure 8 where an intense seismic activity was observed.

The map shown in Figure 12 illustrates the spatial distribution of b-value for northern Algeria. In order to avoid any saturation effect, only 1000 earthquakes extracted from the catalog were used to build this map. The color bar attached to this map emphasizes at the same time the non-homogeneity of this distribution in terms of values of b and geographical location. Two zones corresponding to highest values of b are located in the northeast and in the central part of the northwest regions. Knowing that high values of b indicate an increase of the number of small events, these two zones are therefore highly seismogenic but with weak earthquakes. The high values of b are highlighted by dark blue color. This indicates an intense geological tectonic activity and consequently an accumulation of stress which can induce a likely occurrence of strong earthquakes.

11. b-Value Variations vs. Magnitude

Figure 13 shows the variations of the b-value parameter as a function of the magnitude. The graphic representation exhibits fluctuations of the values of b. However, this fluctuation seems ordered since one can detect a form of linearity with respect to a range of magnitudes.

Figure 11. b-value variations vs. time.

Figure 12. Spatial distribution of b-value.

12. Fractal: Theoretical Background

The concept of fractal analysis has become a widely used tool in the descriptive, studies and understanding of some geophysical problems (Turcotte, 1981) . This theory essentially encompasses so-called fractal geometry and the fractal scaling property. Fractal geometry was initially formulated by the Franco-American mathematician B. Mandelbrot in terms of scale invariance and self-similarity observed in certain objects or phenomena (Mandelbrot, 1982; Crownover, 1995; Feder, 1989) . In addition to these two characteristics, we must add the fundamental parameter of fractal dimension which will play a key role in this study.

Figure 13. b-value variations vs. magnitude.

Fractal geometry provides a tool for investing at different scales of fractal objects and processes. It also allows the observation of the constitutive elements of the irregular shapes which can not be revealed by the Euclidean conventional geometry. The major factor of the fractal theory is undoubtedly the fractal dimension $D_F$ as well as its many variants including the capacity dimension or box counting fractal dimension and the correlation dimension $D_{CR}$ (Klinkenberg, 1994) . A self-similar object is identical or approximately similar to any subset of itself which can be continuously subdivided into ever smaller elements whose structure is a scale-copy of the whole. A seismic regime can be assimilated to a set of earthquake swarms distributed in space and time. The grouping of these swarms gives a representation similar to an aggregate whose size and distribution follow the relationship

$N\left(r>R\right)=C{R}^{-D_F}$ (4)

where $N\left(r>R\right)$ is the full number of objects whose size r is greater than a reference size R and $D_F$ is the fractal dimension. Equation (4) shows that the value of the fractal dimension $D_F$ depends on both the shape of the individual object and the spatial scales that govern the events distribution. By extension, variations in the fractal dimension could be an indicator of the spatial heterogeneity of the distribution. A high value of the fractal dimension implies an increase in the number of objects of small size and hence of the heterogeneity. The fractal dimension thus becomes an indicator of heterogeneity of the spatial distribution of events. Several other approaches to the definition of the fractal dimension have been proposed by various authors.

13. b-Value vs. Fractal Dimension

Seismicity through its different analysis parameters highlights fractal structures in space, time and magnitude distribution as described by the fractal dimension of the epicentre distribution. Among these parameters the b-value, deduced from the Gutenberg-Richter relationship, plays a key role in the analysis of seismological data. Indeed, b-value has the advantage of being correlated with several other parameters of seismicity of which the fractal dimension (Hirata, 1989) . One of these relations was that proposed by Aki (Aki & Richards, 1980) and formulated as follows

$D_F=\frac{3b}{c}$ (5)

where c defines a constant closely related to the size of the seismic events:

・ For small events: $c=1$ so $D_F=3b$ .

・ For intermediate events: $c=1.5$ so $D_F=2b$ .

・ For large events: $c=2$ so $D_F=3b/2$ .

Generally an, with some exceptions, Algeria’s seismicity catalog shows intermediate-size earthquakes. The second case cited above may be accepted i.e. $D=2b$ , the value of b being that provided by the Gutenberg-Richter law as shown in Figure 10: $b=0.773$ therefore $D_F=1.546$ . Smalley et al., 1987 have only retained the second case $D_F=2b$ which gave, in the context of the seismicity map of Algeria, a perfectly acceptable fractal dimension.

14. Fractal Correlation Dimension

There are various methods to compute the fractal dimension of a discrete space or time series apparently randomly distributed. The most popular and the first one used among these methods is the so-called “box counting” or capacity fractal dimension. As introduced by B. Mandelbrot (Mandelbrot, 1982) , it consists in counting the number of small squares (in two dimensions) or cubed (in three dimensions) required to cover the whole of the fractal object under consideration. Therefore, this number is compared with the used box size. The limit of the ratio when the size of a square or cube tends towards 0 defines the box counting dimension noted $D_{CB}$ (Equation (6))

$D_{CB}=\underset{r\to 0}{\lim }\frac{\log \left[N\left(r\right)\right]}{\log \left(\frac{1}{r}\right)}$ (6)

where $N\left(r\right)$ is the number of no, empty boxes (squares or cubes of size r).

The first disadvantage encountered when applying the capacity dimension is that the latter is purely geometric and depends on the density of the number of points of the spatial or temporal series. In the case of a spatial distribution of seismic events, this distribution can be assimilated to a dynamic process. The capacity dimension does not take this dynamic into account, in particular the distance that can separate two samples. To correct this drawback of the box-counting dimension, another empirical dimension called fractal correlation dimension has been proposed. The assessment of this dimension is based on the set of sample points $\left(Z_1,Z_2,\mathrm{...},Z_n\right)$ belonging to the considered fractal object. For a chosen radius r “sufficiently small”, we compute the proportion of pairs $\left(Z_i,Z_j\right)$ denoted $C\left(r\right)$ such as their inter-distance is less than r. This procedure of computing is repeated for many values of r, say $\left(r_1,r_2,\mathrm{....},r_n\right)$ . These values must be chosen sufficiently small. $\ln \left[C\left(r\right)\right]$ is then regressed by $-\ln \left(r\right)$ . The slope of the regression line is the $D_{CR}$ correlation dimension, defined as

$D_{CR}=\underset{r\to 0}{\lim }\left[-\frac{\ln \left(C\left(r\right)\right)}{\ln \left(r\right)}\right]$ (7)

Applied to the data of the seismological catalog this procedure of calculation of the correlation fractal dimension leads to the value 1.4 +/− 0.03 i.e. $1.61<D_{CR}<1.67$ . It was previously calculated the fractal dimension using the b-value. The result was $D_F=1.546$ . The two fractal dimensions are therefore substantially equal (Figure 14).

15. Geostatistical Analysis Based Semivariogram

The statistical analysis by application of the so-called statistical variography method consists mainly of the definition of an experimental variogram calculated from the data and then of the selection of a theoretical variogram model adapted to the data to be processed (Yarus & Chambers, 1994) . The variogram belongs to the class of descriptive statistics techniques (Amorese et al., 2010) . It graphically illustrates and characterizes the spatial continuity of the data set i.e. its roughness. The experimental variogram is calculated by averaging one-half the squared differences of the variable under consideration on all pairs of observations with the specified value separation distance and direction (Equation (8))

$\gamma \left(u,h\right)=\frac{1}{2N\left(a,h\right)}\text{}{\displaystyle \underset{i=1}{\overset{N\left(a,h\right)}{\sum }}{\left[Z\left(i\right)-Z\left(i+h\right)\right]}^2}$ (8)

Figure 14. Estimation of the fractal correlation dimension.

where $\gamma \left(u,h\right)$ is the semi-variogram estimated for the distance h, $Z\left(i\right)$ and $Z\left(i+h\right)$ are the samples of the variable at the points i and $i+h$ , $N\left(u,h\right)$ defines the number of pairs of data observed in the direction $u$ separated each other by distance lag h. When the direction u is not specified, the following model is then obtained

$N\left(h\right)=\frac{1}{2N\left(h\right)}{\displaystyle \underset{i=1}{\overset{N\left(h\right)}{\sum }}{\left[Z\left(i\right)-Z\left(i+h\right)\right]}^2}$ (9)

16. Models and Components of Variograms

According to the data, the theoretical variogram is selected from a set of other several variograms so as to coincide as much as with the experimental variogram (Figure 15).

The seismological data of this study require the choice of an exponential variogram (Figure 16) whose equation is defined as follows

$\gamma \left(h\right)=\{\begin{array}{l}{C}_0+C\left[1-\exp \left(-\frac{3h}{a}\right)\right]\text{if}h\ne 0\\ 0\text{if}h=0\end{array}$ (10)

where $C_o\ge 0,C>0$ and $a>0$ .

The exponential variogram model reaches its sill fairly rapidly and asymptotically. For this model, the range is the value of the distance lag where the variogram reaches 95% of the sill. An example of an exponential variogram model is illustrated by Figure 16 which carries its main descriptive parameters.

A variogram is characterized by at least three descriptive parameters cited as follows: the nugget $C_0$ , the sill identified at the ordinate $C+C_0$ and the range a which represents the lag distance in which the plateau reaches its horizontal shape (Kouadri et al., 2012; Barnes, 1991) .

A non-objective analysis of each of these parameters may the quality of the empirical variogram interpretation. Concretely each parameter defines a specific characteristic of the variogram. The nugget $C_0$ or intercept of the variogram could be attributed to measurement errors or a scale spatial variation (Figure 16). In the case of a continuous variable, the nugget effect should be zero since at

Figure 15. Usual models of semivariograms.

Figure 16. A typical model of exponential variogram.

zero distance and the samples will be the same. However, some variables may change suddenly. A difference over a very short distance is then observable.

When analyzing a semivariogram, it appears at a certain distance, a stabilization of the model. The distance at which the model begins to flatten is called range noted a. If sample locations are separated by a distance smaller than the range, then the samples are spatially correlated. In the opposite case, the samples will be spatially uncorrelated. The value of the semivariogram corresponding to thee range is called sill expressed by $C+C_0$ . Mathematically, the sill is defined

as e $\underset{h\to +\infty }{\lim }\gamma \left(h\right)$ of the random field of data. In another words, the sill is simply variance of the variables.

17. The Variogram as a Geostatitical Analysis Tool

Depending on the statistical objectives and spatial characteristics of the data under consideration a more expressive variogram may be required. For a uniform distribution of the seismic events, a unidirectional or isotropic variogram may be sufficient to describe and to understand that distribution. However in the case of random or irregular distribution, a directional or anisotropy variogram may be used to assess the degree of anisotropy and to extract the major direction of spatial continuity. Anisotropy can be identified by the covariance function or by the nature of the experimental variogram model used for the geostatistical analysis. In practice, a suspected spatial anisotropy can be carried out by the high variability displayed by the experimental variogram which becomes thus disturbed. A directional variogram has the property of varying in structure and fractal dimension by a preliminary choice of a direction of analysis. This behaviour may describe the anisotropic tendency of the data. It also measures the possible dissimilarities expressed as a function of the separation distance on which the variogram is conceived. Indeed, the variable under consideration can change spatially more rapidly in a direction rather than in another.

18. Variogram Approach Based Analysis of Algeria Catalog of Seismicity

Starting from the declustered seismic catalog, the first step is to draw the so-called experimental variogram model. The shape of the graph of the experimental variogram is then compared with a graph of the theoretical variograms in order to extract the mathematical model. As shown in Figure 17, this variogram is represented by a black continuous line. It is obvious that the obtained model is of the exponential type.

The analysis of Figure 17 provides the following numeric values of the exponential semivariogram parameters :

・ Nugget: $C_0=0.13160$ .

・ Sill: $C_0+C=0.42820$ .

・ Range (estimated): $a\approx 0.64$ .

The graph of the exponential semi-variogram describing the spatial distribution of the data of the seismicity catalog of Algeria thus has an intercept represented by the non-zero nugget $C_0$ . This is a consequence of the reliability to be given to the data of the catalog. Indeed, the conversion of seismic data from the intensity scale to the scale of magnitudes can be a source of errors. Let us add to this probably measurement or recording errors.

The calculated sill $C+C_0$ refers to the variance of the set of magnitudes. The variance scale shown in Figure 15 is normalized. Since the variance defines a dispersion parameter of the samples, the value $C_0+C=0.42820$ means that the data contained in the seismic catalog tend rather towards a spatial heterogeneity. The estimated range $a\approx 0.64$ defines the lag distance at which 95% of the variance i.e. 0.4068 is reached. Compared with the value of the sill, it means that beyond this distance the magnitudes are no longer correlable.

Figure 17. Semivariogram of the magnitude data.

19. Semivariogram Cloud Analysis

The variogram cloud is a graph that highlights the average dissimilarity between two observations as a function of the distance separating them spatially. The variogram cloud is applied in the analysis of spatial data to better identify the dissimilarity of the couples and to compensate the deficiencies analysis details not provided by the classical variogram. The variogram cloud plot individualizes the contribution of each data pair to the global variogram (Klinkenberg, 1994; Wackermagel, 2003) . Compared to the classical variogram, it allows to have a wider view of the spatial variation and to identify the unusual behavior of certain points. The latter may be aberrant points that are in fact unnecessary in geostatistical analysis.

The construction of a cloud variogram is based on the calculation of the quadratic difference of the samples relatively to lag distance. The possible dissimilarities are then plotted against separation distance of the pairs of samples. The resulting graph forms the variogram cloud for a fixed class of distance.

20. Algerian Seismicity Catalog vs. Variogram Cloud Analysis

Figure 18 shows a model of variogram cloud calculated for the first class of distance lag. Two distinct areas appear in this figure. A compact zone formed by a gap of pairs of points localized to the weak variances and therefore to strong correlation. The mean curve drawn through this set of points carries out an image similar to that of conventional variogram. The second, sparser area is populated by randomly distributed points, which are probably outlier values that do not participate in the geostatistical study of the spatial distribution of magnitudes

21. Directional Variography

Semivariogram based analysis must not only measure the dissimilarity of the data against distance but also investigate anisotropic data. When using a conventional semivariogram the variable is studied without taking into account the

Figure 18. Semivariogram cloud of the magnitude data.

directions of investigation. In the isotropic variogram model only the modulus of the vector is used (equation). Directional variography introduces another parameter which is the direction of azimuthal investigation (Figure 19). Thus, a study of spatial variability in the distribution of seismological data could be included in geostatistical analysis (Ogata & Katsura, 1993) . It should be noted, however, that in a directional semivariogram only the values of the range and the sill change while the value of the nugget remains constant. In other words, directional analysis can not override measurement errors or those from small scales

22. Anisotropic Variogram Surface

This form of representation of an anisotropic semivariogram also called 2D-variogram map is simply a visual image of semi-variance for various directions. The two coordinate axes contain the separation distances or lags. Thus the abscissa axis corresponds to the E-W direction lags and the Y axis to the N-S direction. The center of the image corresponds to the origin of the semivariogram for each direction. The color bar inserted next to the anisotropic semivariance surface denotes a scale of values for this semivariance. The result of the application of the 2D-variogram map on the data of the seismic catalog of Algeria is illustrated by Figure 20. The direction of the most significant analysis corresponds to $\alpha =45˚$ . Two whitish areas symmetrical with respect to the origin appear clearly. According to the color scale, they indicate high values of the semivariance and therefore low values of the correlation of the seismic events located in this direction. This last corresponds to the geographic direction NE-SW which is precisely the compressive direction of the tectonic regime of the north of Algeria.

23. Variogram vs. Fractals

One of the simplest semivariogram models is the power model given by the following equation

$\gamma \left(h\right)=a{h}^b$ (11)

Figure 19. Anisotropic semivariogram principle.

Figure 20. Surface semivariogram vs separation distances E-W and N-S.

Coefficients a and b are real numbers and they represent respectively a quantification of the rate of variation and a scaling factor. The scaling factor a plays a key role in the shape of the variogram depending on whether $0<a<1$ or $1<a<2$ (Figure 21). These two cases will be applied to an interesting variogram model which is bounded directly to the fractal theory. Its formulation is of the following form

$\gamma \left(h\right)=\frac{1}{2}a{h}^{2H}$ (12)

If we carry out a contraction of scale on the variable h by taking $\frac{h}{k^H}=h^{\prime }$ , the corresponding variogram is given by

$\gamma \left(h^{\prime }\right)=\frac{1}{2}a{\left(\frac{h}{k}\right)}^{2H}$ (13)

Relationship (13) is precisely of the same form as that of the initial variogram. A self-similarity state is thus shown in the case of the variogram. A junction between variogram and fractals is therefore established.

24. Fractal Dimension vs. Fractal Variogram

Two kinds of geometric dimensions can characterize the spatial distribution of data such as those provided by a seismic catalog. These are the natural or topological dimension $D_T$ and the fractal dimension $D_F$ . Generally, the fractal dimension is comprised between two consecutive topological dimensions, it is then fractional. Fractal objects are characterized by their apparent complexity,

Figure 21. First simplest model of semivariogram.

irregularity and random fluctuations. The “fractality” of an object is an intuitive perception of its irregularity and roughness assessed by the fractal dimension. The more an object is irregular, the higher its fractal dimensions. The relationship between the topological dimension $D_T$ and the fractal dimension $D_F$ is provided by the physical process called the Brownian Fractional Movement (fBm), the study of which is outside the scope of this study. Equation (13) is a model of fBm indexed by the parameter H such that $0<H<1$ . The three quantities H, $D_T$ and $D_F$ are unified by the following relationship

$D_F=D_T+1-H$ (14)

For $D_T=1$ we obtain

$H=2-D_F$ (15)

Equation (14) is the basis of the geostatistical analysis by fractal variogram of the seismicity catalog data. Equation (13) can be written as

$\gamma \left(h\right)=\frac{1}{2}a{h}^{2\left(2-D_F\right)}$ (16)

Bi-logarithmic coordinates applied on Equation (16) i.e. $\log \left[\gamma \left(h\right)\right]$ against $\log \left(h\right)$ provides the equation of a straight line pe $m=2\left(2-D_F\right)$ . The fractal dimension $D_F$ is then easily calculable:

$D_F=2-\frac{m}{2}$ (17)

This procedure will allow the analysis of the spatial anisotropy of the data of the seismic catalog jointly by using the directional semivariogram and the fractal dimension.

25. Spatial Anisotropy Using Directional Variogram and Fractal Dimension

Four a directions of analysis of the spatial anisotropy were chosen. These directions correspond to those most commonly used. By applying the previous calculation procedure, the fractal dimensions btained are summarized by Table 1.

Figure 22 illustrates the set of the directional variograms in bilogarithmic coordinates as well as the lines whose slopes will provide the fractal dimensions according to Equation (15). The calculated fractal dimensions are relatively high, close to the surface topology dimension i.e. 2, whatever the direction considered. There is therefore no preferred anisotropic direction in the distribution of earthquakes. This corroborates the thesis of high seismicity all along the north of Algeria. Let us notice especially the value of the fractal dimension $D_F=2$ which is exactly the value of the topological dimension obtained for the direction $\alpha =90˚$ i.e. for a transverse path.

26. Result and Discussion

The most popular scaling power law for earthquakes i.e. the Gutenberg-Richter’s relation is the basis of this study. It describes for a given seismogenic area the distribution of earthquake occurrence. Two important seismological parameters are provided by this law namely the a-value and the b-value the a-value and especially the b-value which depends on the tectonic regime induced by the stress field. The Gutenberg-Richter’s b-value was the first seismological parameter calculated and interpreted.

Table 1. Fractal dimension vs. azimuthal analysis.

Figure 22. Estimation of the fractal dimension vs. directional anisotropic analysis.

The b-value is of great importance, because the concept of geometrical self-similarity it has been shown that b-value can be directly related to the fractal dimension. The graphically estimated value of the b-value for the seismicity of Algeria is $b=0.773\pm 0.06$ . The analysis of this global value as well as the distribution map of this seismic parameter reveals a moderate, frequent and uniformly distributed seismic activity throughout the north of Algeria.

The concept of fractal and related concepts of self-similarity and scale invariance have been introduced to understand the spatial organization of cluster epicentres that obeys an underlying order known as fractal structure.

The value of the fractal dimension is considered as an index of the complexity that can prevail in a time series such as a data of magnitude. The moderate seismicity of northern of Algeria has suggested using the relationship $D_F=2b$ i.e. $D_F=1.546$ . This value of the fractal dimension is intermediate between the two topological dimensions 1 and 2, first confirming average data complexity and non-uniform epicentre distribution. However, the fractal dimension $D_F$ does not take into consideration the clustering mode in the spatial distribution of epicentres. Another type of fractal dimension is then proposed: the correlation dimension denoted $D_{CR}$ . Graphically the value $D_{CR}=1.64\pm 0.03$ was obtained. Fractal dimension $D_F$ and correlation dimension $D_{CR}$ are substantially equivalent. They are substantially equivalent. The clustering mode of earthquakes in northern Algeria thus presents a not very complex distribution.

Basically, the mean and variance parameters are of great importance in the geostatistical analysis of the time series because they are respectively indicators of the homogeneous and isotropic character of the data. With respect to homogeneity, the variographic approach was used in this study because it provided an experimental variogram of exponential type with a non-zero nugget $C_0$ which gives some reliability to the data despite the likely uncertainties of measurement or recording. The calculated sill $C+C_0$ means that the data contained in the seismic catalog tend rather towards a spatial heterogeneity.

The study of the spatial anisotropy of the data was carried out by means of the variogram which showed an NW direction of anisotropy i.e. the same direction of the compressive geodynamic regime previously discussed.

The junction between semivariogram and fractal dimension is established through the directional anisotropic variogram. For a chosen range of directions, the fractal dimension was calculated; thus for the direction $\alpha =45˚$ it was obtained the fractal dimension $D_F=1.942$ close to the topological dimension $D_T=2$ . This means that the density of the earthquakes is very high in this direction, which therefore has a high seismic activity.

27. Conclusion

Despite the constraints imposed by the nature and quality of the seismological data, the catalog of seismicity of Algeria was proved to be a fruitful source of information on the seismological activity. The extraction and processing of this information requires the development of an analysis protocol. Statistical analysis in general and the application of geostatistical methods in particular are tools adapted to these objectives. Despite their relevance, current seismological parameters such as b-value can not suffice to carry out an extended survey of seismological events. They lack spatial dimension that is to say a description of the distribution of epicentres. The concepts of conventional or isotropic variogram and its variants such as the variogram cloud, the variogram surface and the directional variogram are additional analytical tools recommended for catalog data analysis. Fractal theory and fractal geostatistics are among the most recent orientations for the analysis of seismological data in the catalog. The first one provided the concept of fractal variogram across the fractal dimension which thus becomes a seismic parameter suitable for the study of the anisotropy of spatial distribution of epicentres for the study area.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Aitouche, M. A., Djeddi, M., & Baddari, K. (2012). Fracal Variogram Based Aftershock Sequences Analysis-Case Study: The May 21, 2003 Boumerdes Algeria Earthquake MW = 6.8. Arabian Journal of Geosciences, 6.
[2]	Aki, K., & Richards, P. (1980). Quantitative Seismology Theory and Methods. San Francisco: H. Freeman Co.
[3]	Amorese, F. W. H. D., Grass, J.-R., & Rydelek, P. (2010). On Varying b-Value with Depth: Results from Computer Intensive Tests for Southern California. Geophysical Journal International, 180, 347-360.
[4]	Barnes, R. (1991). The Variogram Sill and the Sample Variance. Mathematical Geology, 23, 673-678. https://doi.org/10.1007/BF02065813
[5]	Bellalem, F., Mobarki, M., & Talbi, A. (2008). Analysis of Spatio-Temporal Seismic Activity in in Northern Algeria. Proceeding of the 14 World Conference on Earthquake Engineering. Beijing.
[6]	Benouar, D. (1994). Materials for the Investigation of the Seismicity of Algeria and Adjacent Regions during the Twentieth Century. Annals of Geophysics, 37, 459-860.
[7]	Carminati, C., Doglioni, C., & lastrino, M. (2012). Geodynamic Evolution of the Central and Western Mediterranean Tectonic vs. Igneous Petrology Constraints. Tectonophysics, 579, 173-192. https://doi.org/10.1016/j.tecto.2012.01.026
[8]	Cressie, N. A. (1993). Statistics for Spatial Data. New York: Wiley & Sons.
[9]	Crownover, R. (1995). Introduction to Fractals and Chaos. Burlington: Johns and Barlett Publishers.
[10]	Feder, J. (1989). Fractals. New York: Plenum Press.
[11]	Goovaerts, P. (1997). Geostatistics for Natural Resources Evaluation. Oxford: Oxford University Press.
[12]	Gutenberg, B., & Richter, C. F. (1942). Earthquake Magnitude Intensity and Acceleration. Bulletin of the Seismological Society of America, 32, 163-191.
[13]	Hirata, T. (1989). A Correlation between b Value and Fractal Dimension of Earthquakes, JGR Solid Earth, 94, 7507-7514. https://doi.org/10.1029/JB094iB06p07507
[14]	Kenzi, Mc. D. P. (1972). Active Tectonic of the Mediterranean Region (Soc. 30). Geophy. J.R. Astronomy.
[15]	Klinkenberg, B. (1994). A Review of Methods Used to Determine the Fractal Dimension of Linear Features. Mathematical Geology, 26, 23-46.
[16]	Kouadri, A., Aitouche, M. A., & Zelmat, M., (2012). Variogram Based Fault Diagnosis in an Interconnected Tank System. ISA Transactions.
[17]	Legrand, D. (2002). Fractal Dimensions of Small, Intermediate and Large Earthquakes. Bulletin of the Seismological Society of America, 92, 3318. https://doi.org/10.1785/0120020025
[18]	Mandelbrot, B. B. (1982). The Fractal Geometry of Nature (2nd ed.). Londra: W. H. Freeman.
[19]	Mathéron, G. (1970). La théorie des variables régionalisées et ses applications, Cahiers du Centre de Morphologie Mathématique de Fontainebleau. Paris: Ecole des mines.
[20]	Mc Nutt, S. R., Sanchez, J. J., & Wyss, M. (2004). Spatial Variations in the Frequency-Magnitude Distribution of Earthquakes at Mount Pinatubo Volcano. Bulletin of the Seismological Society of America, 94, 430-438.
[21]	Ogata, Y., & Katsura, K. (1993). Analysis of Temporal and Spatial Heterogeneity of Magnitude-Frequency Distribution Inferred from Earthquake Catalogues. Geophysical Journal international, 113, 727-738. https://doi.org/10.1111/j.1365-246X.1993.tb04663.x
[22]	Olea, R. A. (2009). A Practical Prime on Geostatistics. USG Geology Survey.
[23]	Schaefer, A. M., Danielle, J. E. & Wenzel, E. (2014). Application of Geostatistical Methods and Machine Learning for Spatio-Temporal Earthquake Cluster Analysis. American Geophysical Union Fall Meeting.
[24]	Smalley, R. F., Chatelain, J.-L., & Prevot, R. (1987). A Fractal Approach to the Clustering of Earthquakes. Bulletin of the Seismological Society of America, 62.
[25]	Turcotte, D. L. (1981). Fractals and Chaos in Geology and Geophysics. Cambridge: Cambridge University Press.
[26]	Wackermagel, H. (2003). Variogram Cloud. In Multivariate Geostatistics. Berlin: Springer.
[27]	Wierner, S. (2001). A Software Package to Analyse Seismicity ZMAP. Seismological Research Letters, 72, 373-382.
[28]	Wiemer, S., & Wiss, M. (2002). Mapping Spatial Variability of the Frequency-Magnitude Distribution of Earthquakes. Advances in Geophysics, 45, 1-40.
[29]	Wierner, S. (1996). ZMAP User Guide. XYZ Corporation.
[30]	Yarus, S., & Chambers, M. A. (1994). Stochastic Modelling and Geostatistics: Principles, Methods and Case Studies (No. 3). AAPG Computer Applications in Geology.