ABSTRACT

Probabilistic seismic hazard assessment (PSHA) takes into account as much data as possible for defining the initial seismic source zone model. In response to this, an algorithm has been developed for integration of geological, geophysical and seismological data through a spatial index showing the presence or absence of a potential seismic source feature in the input data. The spatial matching index (SMI) is calculated to define the coincidence of independent data showing any indications for existence of a fault structure. It is applied for hazard assessment of Bulgaria through quantification of the seismic potential of 416 square blocks, 20 × 20 km in size covering the entire territory of Bulgaria and extended by 20 km outside of the country borders. All operations are carried out in GIS environment using its capabilities to work with different types of georeferenced spatial data. Results show that the highest seismic potential (largest SMI) is observed in 56 block elements (13% of the territory) clearly delineating cores of the source zones. Partial match is registered in 98 block elements when one of the features is missing. Not any evidence for earthquake occurrence is predicted by our calculation in 117 elements, comprising 28% of the examined area. The quantitative parameter for spatial data integration which is obtained in the present research may be used to analyze information regardless of its type and purpose.

Keywords:

Seismic Hazard, Data Integration, Matching Index, Spatial Analysis, Bulgaria

1. Introduction

Probabilistic Seismic Hazard Assessment (PSHA) requires compilation of a knowledge base which is formed by the identification and documentation of relevant possible earthquake sources in a region. A key part of this analysis is the identification of regional or local seismotectonic domains that, based on geologic, geophysical, and seismological information, may be interpreted to have relatively consistent spatial and temporal variations in historical seismicity.

The territory of Bulgaria can be assigned to one of the most seismically active and rapidly deforming regions within the continents. As a result of the neotectonic evolution and dynamics of the Aegean extensional system a number of seismological structures exist and accommodate the crustal deformation in Bulgaria [1] thus posing serious challenges to seismic hazard assessment. Previous studies relied mainly on a visual correlation of the lineaments obtained from different types of data. Reference [2] presents a complex analysis of the existing seismological, geological, geophysical and geodetic information using the results obtained by a big team of experts (more than 15 persons) who independently traced the lines on the separate maps and after a visual inspection determined the position of the complex lineaments. Afterwards, a successive generalization was performed and according to predefined rules the outlined seismogenic lineaments were ranked to a certain level of confidence and accepted for the final map compilation. A huge amount of time, people and resources have been invested to ensure objective and reliable integration of the input data.

For the seismic hazard assessment accomplished in 2009, the seismotectonic map was elaborated [3] through digital juxtaposition of geological, geophysical and seismological maps using direct visual evaluation of their coincidence. Although this technique allowed more precise comparison of the fault indicators, it remained still possible to underestimate or overestimate some evidences important for the accuracy of the seismotectonic model.

To overcome the weaknesses of the above described methods we propose a newly developed index which shows the existence or not of a match in the spatial distribution of any two features (characteristics) regardless their type or units of measure. For the purposes of this study we compare the spatial distribution of active faults, Bouguer gravity lineaments derived from the gravity anomalous field and seismological data by means of earthquake’s space-temporal parameters. We process and integrate the available data sets in GIS and asses the seismogenic potential of the Bulgarian territory calculating a specially developed matching index thus avoiding subjective evaluation of a visual comparison and the danger of omissions.

2. Data

During the last two decades a number of high quality geological, geophysical and seismological data were collected and used for the purposes of seismic hazard assessment. To perform quantitative spatial analysis supporting the definition of the studied region seismotectonic characteristics we compile an input database using the following three independent geo data sets: geological, geophysical and seismological. The territory of Bulgaria is divided into square blocks of size 20 × 20 km. It is extended by 20 km outside of the country borders to allow coverage of data from the nearest border region that have an impact on the seismic hazard as well (Figure 1).

2.1. Geology

The first key dataset which is used in the present study is the spatial distribution of active faults (AF) on the territory of Bulgaria. Data (Figure 2(a)) are modified from the map of active faults that is presented in [3].

The faults are classified into three categories depending on their characteristics and degree of knowledge: active, potentially active and possibly active (with unverified activity). In the category of active faults are those for which there are clear evidences of the Late Pleistocene or Holocene activity. To the category of potentially active faults are referred those Faults marked on the base of indirect information fragments or lines that can be activated in the contemporary tectonic regime. As possibly active, with unproven activity are identified faults or lineaments with scarce data on their ruptured origin and Late Pleistocene or Holocene activity.

All classes of the active faults described above are included in the study because we use a conservative hazard assessment approach in our analysis.

The length and end-points’ coordinates of the outlined faults are specified in a data set table (Figure 2(a)). More than 130 structures with length larger than 5 km are included. The length of the longest fault is 68 km (near Plovdiv), while the average length of structures is 18 km. The dominant fault strike is in NW-SE direction.

2.2. Geophysics

A variety of geophysical methods and data could be used to provide evidences about geological structures in depth. Most appropriate for regional interpretation are geomagnetic and gravimetric anomalies, which allow application of specialized transformations to increase the direct correlation between the observed anomalous field and the respective geological source.

As more suitable for detection of fault-like features we chose the modulus of the total horizontal gravity gradient (THG) due to his ability to highlight transition type anomalies (gravity lineaments―GL) which are connected with faults, horst or graben structures, block borders, etc. [4]. THG is calculated from the Bouguer anomalous field map of the Bulgarian territory [5].

The anomalies’ intensity is in the range of 176 mGal from negative values in the Rila-Rhodope area to positive values at the Black Sea coast. The average gravity gradient from west to east is 0.32 mGal/km. The calculated local total horizontal gradients show values up to 8 mGal/km.

The horizontal derivatives along two orthogonal axes are calculated and geometrically summed in a grid with a 3 km. cell size over the entire territory of Bulgaria [6]. The map of THG lineaments indicates the axes of pronounced gravity transitions by the lines of maximum gradient values. The most intensive among them are marked by long black lines (Figure 2(b)). The length and coordinates of the edges of lineaments are specified in a data set table. The number of marked gravity transitions is more than 110. The longest transition is 42 km of size (to the northeast of Sofia), and the average length of the transitions found is about 19 km. The dominant orientation of their strikes is in WNW-ESE direction following the main lineaments of the known structures of the lithosphere in Bulgaria.

2.3. Seismology

The last but not the least, an earthquake catalogue containing 755 shallow earthquakes with a magnitude M_W > 3.0 occurred in the territory under consideration (Figure 2(c)) is used as the third dataset. Earthquakes cover the period from the 1st century BC up to year 2016 and are based on available historical documents and information from the sources described in [7].

To insure homogeneity in energy domain and compatibility of seismological information magnitude estimates are converted in [8] to the most widely used (in recent years) scale-seismic moment magnitude M_W [9].

The de-clustered and homogenized catalogue was also studied for the completeness (e.g. [10] ) using the Stepp’s test [11] , which is essential for the reliable evaluation of the seismic statistical parameters used in probabilistic seismic hazard assessment.

The described information is organized in a dataset (Figure 2(c)) where each earthquake is specified with space-temporal parameters (date, time, longitude, latitude, depth) and magnitude (D, T₀, φ, λ, h, Mw,). The dataset contains 755 events occurred in and near Bulgaria with the largest magnitude M_W = 7.6 (the 1904 earthquake occurred near Krupnik (23.2˚E, 41.8˚N)).

The parameter used in the present analysis is the logarithm of the seismic moment M₀ calculated from the moment magnitude scale [12] using the equation:

$M_w=\frac{2}{3}\log M_0-10.7$ (1)

were M_w is the moment magnitude from the seismological dataset described above.

The seismic moment M₀ is a scalar value depending on the area of fault rupture, the average amount of slip (displacement), and the shear modulus of the rocks involved in the earthquake. M₀, [Nm] is in the basis of the moment-magnitude scale introduced by Kanamori [13] which is often used to compare the size of earthquakes.

3. Method

The considered territory is gridded with a step of 20 km. As a result 416 square blocks are obtained each of them containing geological, geophysical and seismological data indicating the presence or absence of a potential seismogenic features. We analyze available information in pairs (geological and geophysical; geological and seismological; geophysical and seismological) by using a specially derived parameter named “spatial matching index” (SMI).

The technique which we apply in this study consists in calculating of the SMI which assesses the existence or absence of any two V_i and V_j features representing the available information in the input data of each element of the spatial square blocks grid. The SMI, marked by the symbol Q in the Equation (2), for a given element with a cell number n is estimated using the formula:

$Q_n\left(V_i,V_j\right)={\left[\frac{V_i+V_j}{\sqrt{V_i^2+V_j^2+1}}\right]}_n,n=1,2,3,\cdots ,k,\text{and}i\ne j$ (2)

where variables V_i and V_j might be different geological, geophysical or seismological characteristics having the same or different measurement units. They are derived from the available information and prescribed as non-negative numbers for each square element (k = 416) of the used grid.

According to the obtained value of Q is assessed the existence or not of a match in any of the three feature pairs (characteristics) V_i and V_j. The limits of the matching index Q are strictly defined and used in the analysis as follows:

- if V_i > 0 and V_j > 0, then Q is within the limits [1.0; $\sqrt{2}$ ), and the grid element with number i; is marked by “existence” of both characteristics;

- if V_i = 0 and V_j = 0, then Q = 0, the grid element is marked by a “zero value matching” or “missing” for both characteristics;

- if V_i > 0 and V_j = 0, or V_i = 0 and V_j > 0, i.e. only one of the two features has an expression larger than 0 in the network element, then the Q value falls within the range (0; 1), with the assessment of “non-match” in the manifestation of the two characteristics.

As is seen in Figure 3, the maximum values of Q are obtained when V_i = V_j, which means that as close are the compared parameters as high the matching index will be.

Above methodology is applied using ArcGIS on the three key datasets: active faults, lineaments of the total horizontal gravity gradient and catalogue of earthquakes (epicenter’s distribution). For each element of the grid were calculated the sum of the active faults’ length ( ${\displaystyle \sum AF}$ ), sum of the gravity lineaments ( ${\displaystyle \sum GL}$ ), and logarithm of the aggregated seismic moment ( $\lg {\displaystyle \sum M_0}$ ). These resulted in 122 single square blocks of the geological data (AF), 157 blocks of geophysical data (GL) and 141 blocks of seismological data (M₀) having V > 0.

After the three mutual matching indices Q_AF-GL (from the geology―geophysics pair), Q_AF-lgM (from the geology―seismology pair) and Q_GL-lgM (from the geophysics―seismology pair)were calculated using Equation (2), it was possible to determine whether or not the three data types are matched in any element of the computing network using an “integrated” spatial matching index Sq:

$Sq=Q_{AF-GL}+Q_{AF-lgM}+Q_{GL-lgM}$ (3)

Four different cases are observed for Sq according to the individual Q values (Equation (2) and Equation (3)):

Sq = 0 “None”―not any feature exists in the cell;

0 < Sq < 2.0 “Non-match”―only one feature exists;

2.41 < Sq < 3.41 “Partial match”―two of the three features exists;

3.0 < Sq < 4.24 “Full match”―all three features are presented in the cell.

4. Results

Results about the juxtaposition of the three independent data sets (geological, geophysical and seismological) obtained from the calculation of an integrated spatial matching index Sq are shown in Figure 4. With the lightest color are marked areas without any indications of presence of a seismogenic structure while the darkest one correspond to a full match of symptoms for a possible earthquake source existence.

The case of non-match means that information for structure existence is coming from only one of the three data sets. Most of these cells are observed in the northwestern Bulgaria due to the presence of gravity lineaments reflating the block borders of the Moesian platform and its coupling with the West Balkan [14]. The lack of earthquakes and observed faults confirms the present tectonic stability of those structures today. Similar is the case southwest of the Burgas region. Due to the number of effusive and intrusive magmatic bodies embedded in the upper crustal section [15] , a lot of corresponding gravity lineaments are delineated while there is not information about earthquakes and active faults.

Partly matching elements of the grid represents areas where two of the three indicators are observed. As such, very clearly stand out the region north-east of Plovdiv (between Stara Zagora and Yambol) and the sub-parallel strip north-west and south of Burgas.

The first one is connected with the densely located fault structures and a lot of earthquake epicenters (see Figure 2), while the second one is due to the presence of faults and gravity lineaments marking the large plutonic bodies with high density and ultra-basic composition located there [16].

In the north-eastern part of the territory the mosaic pattern of the calculated

index is related to large transposed blocks inside the North-Bulgarian arc and neighboring depressions. Block borders are mainly of fault type and a number of epicenters are also spread around.

The most impressing in Figure 4 are the Upper Thrace basin region (around 25˚E, 42˚N) and the entire territory of southwestern Bulgaria. The zone near Plovdiv contains the most intensive gravity gradient anomalies, active faults and well grouped earthquake epicenters. One orientation of the lineaments dominates along Maritsa River outlining the deep Maritsa dislocations, as well as the set of faults in the heterogeneous basement depression. Maritsa deep fault has been active during Upper Cretaceous stage and it became a transmitter of intensive magmatic activity through its multiple fragments in the upper part of the crust. They are part of the Upper Thrace basin which is a neotectonic structure formed between Rhodope massif and Central Srednogorie unit where a cyclic accumulation of Neogene-Quaternery sediments was observed [17]. The present seismotectonic activity is related to these fault systems that generated a large number of earthquakes over the centuries (see Figure 2(c)).

Complicated border faults of the West and Central Srednogorie units (zones according to [14] ), neotectonic structures of the Sofia complex graben and the fault zones along the Struma and Mesta Rivers form the region with the highest observed matching index. Gravity lineaments are caused mainly by plutons and basic intrusions fed from that magma conducting fractures. Geological dataset also contains a great number of dislocations marked as active or potentially active faults (Figure 2(a)) which combined with the high epicenter density result in dark-colored elements in Figure 4.

5. Conclusions

The approach of using GIS for integrating different types of geo data and calculating the spatial matching index (SMI) that we applied in this study clearly outlines the areas with the highest earthquake capacity thus subjective evaluation of a visual comparison, omissions and errors are avoided.

The highest seismic potential (largest SMI) is observed in 56 square block elements which comprise 22,400 km² (13% of the examined territory in and near Bulgaria). The regions near to the cities of Varna, Dulovo, V. Tarnovo, Sofia, Plovdiv, Krupnik possess a full set of studied characteristics. Together with the regions resulted in a partial match of the studied seismic indicators (98 block elements) they form the initial idea of the seismic sources delineation. Lack of evidences for earthquake occurrence is predicted by our calculation for about 28% of the considered area. This result was compared with the seismic hazard map of Bulgaria presented in [18] which is derived from source zones of normal depth for a 475 years recurrence period. The obtained “safe” area corresponds to regions colored in blue, which represents the intensities between 6 and 7 in MSK-1964 scale. Undoubtedly these are values which could be observed as a result from events realized from the neighboring source zones only.

The result from our study is the first attempt to obtain a quantitative parameter for integration of input geo data sets which favor the elaboration of seismic source zones model on the territory of Bulgaria. It might be extended further to cover more information and to result in more accurate evaluations regarding delineation of the seismotectonic domains and seismic hazard assessment.

Acknowledgements

We thank two anonymous reviewers for their critical comments that helped us to improve the manuscript.

The presented work is supported by the WORLD FEDERATION OF SCIENTISTS under the National Scholarship Programme. The financial support for article processing charges is provided by Contract No. D01-161/28.08.2018 (Project “National Geoinformation Center (NGIC)”).

Conflicts of Interest

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

Cite this paper

Trifonova, P., Metodiev, M., Stavrev, P., Simeonova, S. and Solakov, D. (2019) Integration of Geological, Geophysical and Seismological Data for Seismic Hazard Assessment Using Spatial Matching Index. Journal of Geographic Information System, 11, 185-195. https://doi.org/10.4236/jgis.2019.112013

References