Probabilistic seismic hazard in terms of intensities for Bulgaria and Romania – updated hazard maps

G. Leydecker* , H. Busche* , K.-P. Bonjer2, T. Schmitt1, D. Kaiser1, S. Simeonova3, D. Solakov3, and L. Ardeleanu4 1Federal Institute for Geosciences and Natural Resources, Hannover, Germany 2Geophysical Institute University of Karlsruhe, Germany 3Bulgarian Academy of Sciences, Geophys. Instit., Seismolog. Depart., Sofia, Bulgaria 4National Institute for Earth Physics, Bucharest, Romania * formerly at: Federal Institute for Geosciences and Natural Resources, Hannover, Germany


Introduction
Romania and Bulgaria are situated in the Carpathian-Balkan Region as a part of the Alpine-Himalayan seismic belt.They are characterized by high seismicity, and are exposed to a high seismic risk.Over the centuries, the two countries experienced strong earthquakes.Because both countries are now members of the European Union, they are committed to introduce the EUROCODE 8 (EC8).For the application of this building code seismic hazard maps with a single hazard parameter such as macroseismic intensity, peak ground acceleration or spectral acceleration, each is possible, have to be calculated.Based on these results, hazard zones with dedicated referenced peak ground accelerations have to be determined in a further process.In this paper we present a special way to compute such hazard maps.We prefer to use here instead of magnitude the macroseismic intensity to calculate hazard maps, because this parameter is the original parameter Published by Copernicus Publications on behalf of the European Geosciences Union.

M2
Fig. 1.Earthquake epicenter map (I 0 = epicentral intensity) with the seismic source zones.Codes of source zones and their parameters are listed in Table 1.
in historical earthquake catalogues.A great advantage of using intensities is that the very irregular pattern of the attenuation field of the Vrancea intermediate depth earthquakes can be estimated from detailed macroseismic observations.The seismic hazard is computed with EQRISK (McGuire, 1976), based on the probabilistic approach developed by Cornell (1968).
During the last decades, several papers have been published about the seismic hazard in Bulgaria and in Romania.An overview about these approaches is already given in Ardeleanu et al. (2005) for Romania and in Simeonova et al. (2006) for Bulgaria.The presented earthquake hazard maps are an update respectively improvement of the maps in Ardeleanu et al. (2007).

Earthquake catalogue and seismic source zones
The main base of our probabilistic analysis is the earthquake catalogue for SE-Europe (Shebalin et al., 1998) in combination with the Romanian catalogues (for references see Oncescu et al., 1999) and some further national and regional catalogues (Christoskov et al., 1979;Grigorova et al., 1979;Solakov and Simeonova, 1993;Papazachos et al., 2000Papazachos et al., , 2005)).From these different catalogues Ardeleanu et al. (2005) and Simeonova et al. (2006) compiled a homogeneous earthquake catalogue.
As a Poisson model for the seismicity is assumed, only independent events should be included in the analyses.The fore-and aftershocks were removed by using a space-time window.In the case of more than one earthquake within 10 days and 50 km distance only the strongest event is used for the statistics; the others are considered as statistically dependent and rejected from the catalogue.The resulting earthquake catalogue of the considered area consists of nearly 2750 events with epicentral intensities I epic ≥4.
The completeness of the catalogue was examined using Stepp's (1971) test.Stepp's (1971) analysis provides a method for the analytical determination of time intervals in which a particular intensity class (or magnitude) is likely to be reported completely.The test results imply that it is possible to create homogeneous data samples by determining intervals over which earthquakes in different intensity classes are completely reported.
The spatial pattern of seismicity in an area of about 200 km around the territory of Romania and Bulgaria is shown in Fig. 1.The figure shows the earthquakes with I epic ≥4 reported in the time 342 BC-1990 AD, together with the seismic source zones.Each zone is characterized and defined by its own specific seismicity, and geological and tectonic development.From our seismotectonic analysis of the considered area this seems more appropriate than to use specific linear fault structures.It is assumed that an earthquake can occur within a source zone at every point with equal probability.The boundaries of these zones were chosen to reflect the seismicity adequately regarding tectonic units and lithospheric structure.The seismic sources are mostly the same as in Simeonova et al. (2006), and in Ardeleanu et al. (2007), with small differences in Bulgaria (Report GPh.I., 2007) and for Vrancea intermediate depth seismic zone (see below).All earthquakes outside the designed sources are assigned to the "background seismicity".
From the analysis of the depth distribution it was recognized that the earthquakes in all zones occurred in the Earth's crust (h<60 km), with the exception of the events in the Vrancea (Romania) intermediate depth zone.The seismogenic depth of each zone is defined as the depth of maximum seismic energy release and is given in Table 1.
The seismicity within the Vrancea (Romania) region consists of two depths domains: normal depth (less than 60 km) and intermediate depth (60-180 km) earthquakes.Due to extreme irregularities of the isoseismals of intermediate depth earthquakes their effects on seismic hazard have to be treated separately.The size of the Vrancea intermediate depth seismic zone is confined to an area of 20×60 km 2 (Bonjer et al., 2005(Bonjer et al., , 2008)).For each seismic source the intensity-frequency relation was calculated and a maximum possible earthquake was fixed.We assume a truncated exponential distribution as a seismicity model for each seismic source zone, defined by the parameters a and b of the intensity-frequency relation and by the maximum possible earthquake with intensity I max .The parameters a and b of each seismic source zone Nat.Hazards Earth Syst.Sci., 8, 1431Sci., 8, -1439Sci., 8, , 2008 www.nat-hazards-earth-syst-sci.net/8/1431/2008/End year for all regions is 1990; usually, only events with intensities of 5.5 or more are considered.
b) The regression curve in the cumulative intensity-frequency relation for region AF is calculated without the single event with intensity 10.0 MSK for statistical reasons; the next lower class IX contains 9 events.c) For the two regions BA and DA, the events are added and a common intensity-frequency statistic is done, resulting in a common b-value.The a-value is computed using the cumulative number of events with intensity class 7.0 MSK and more for region BA; and with intensity class 6.0 MSK and more for region DA.d) For the two regions BD and PD, the events are added and a common intensity-frequency statistic is done, resulting in a common b-value.The a-value for each region is computed, using the cumulative number of events with intensity class 5.0 and greater.
e) The regression curve in the cumulative intensity-frequency relation for region D2 is calculated without the single event with intensity 10.0 MSK, for statistical reasons; the next lower class IX contains 17 events.
f) The regression curve in the cumulative intensity-frequency relation for region GO is calculated without the single event with intensity 10.0 MSK, for statistical reasons; the next strongest event has intensity 8.0 MSK.g) Only events with intensities class 7.0 or more are considered.
h) The regression curve in the cumulative intensity-frequency relation for region KR is calculated without the two events with intensity 10.0 MSK, for statistical reasons; the next strongest event has intensity 8.0 MSK.i) For the two regions N1 and N2, the events are added and a common intensity-frequency statistic is done, resulting in a common b-value.However, the a-value for each region is computed separately, using the cumulative number of events with intensity class 6.0 and greater.j) The regression curve in the cumulative intensity-frequency relation for region SH is calculated without the single event with intensity 10.0 MSK, for statistical reasons; then the strongest event has an intensity of 8.0 MSK.The statistic was done for a bigger area (thin lines in Fig. 1) because of suspected mislocated events; the used size of the source region SH is shown by thick lines.k) For the regions M1 and M2 as well as for the regions SL and SU, the events are added and a common intensity-frequency statistic is done, resulting in a common b-value for each pair.The a-value for each region is computed, using the cumulative number of events with intensity class 6.0 and greater.
l) The regression curve in the cumulative intensity-frequency relation for region SW is calculated without the single event with intensity 10.0 MSK, for statistical reasons; the next strongest event has intensity 8.0 MSK. are calculated by the least square method using the following equation: where N (I ) is the cumulative number of earthquakes with intensity I ≥I epic .Half intensity values were used, in contrast to the earlier papers (Ardeleanu et al., 2005(Ardeleanu et al., , 2007;;Simeonova et al., 2006) where the half intensities were rounded up to the next full number.The beginning of the time period is chosen according to an assumed completeness of the catalogue for events greater than the lowest intensity values used in the statistics (Table 1).Some of the seismic zones (as for example the seismic sources N1 and N2 in Northern Bulgaria, see Table 1) display low and dispersed seismicity.Since it is not possible to determine a reliable intensity-frequency relation for each of these sources separately, they are treated in a special way (see i) in Table 1) to calculate the intensity-frequency statistics.
All model parameters for PSHA are summarized in Table 1.

Shallow earthquakes
We implemented into the program EQRISK the intensity attenuation function of Sponheuer (1960) which is based on  (1978).The yellow star represents the rupture initiation (epicenter), the white star the rupture termination (stopping) according to Müller et al. (1978).Kövesligethy (1907).The attenuation is given by the following equation: I epic is the epicentral intensity, I site the intensity at a site at hypocentral distance r (km); h is the focal depth (km), and a the absorption coefficient.Here we considered a=0.002 km −1 .For a hazard curve at a site EQRISK cuts all source regions into finite ring segments with their assigned statistical parameters.Then the site intensities caused by earthquakes of each segment are calculated assuming a standard deviation of half an intensity for the attenuation function.The sum of the contributions of all regions finally leads to the annual probability of exceedance at each grid point of the hazard map.These calculations are done for all points between 39 • N to 49 • N and 19.5 • E to 30.0 • E every 0.1 • in latitude and 0.2 • in longitude.Figure 5 shows the seismic hazard for Romania and Bulgaria due to the shallow earthquakes for a recurrence period of 475 years.

Vrancea intermediate depth earthquakes
The macroseismic field of intermediate depth earthquakes of the Vrancea zone Vi (Romania) is very irregular as shown in Figs. 2 and  12.

18.
12. 12. 1977.These figures demonstrate as well as the macroseismic field of the strong earthquake in 1986 (Radu et al., 1987) that the attenuation in the North-West and South-East directions is substantially higher than in the North-East and South-West directions.To account for the apparent azimuthal dependency of the Vrancea intermediate depth earthquake attenuation, different approaches were made and several attenuation relationships were developed (e.g.Mârza and Pântea, 1995;Lungu et al., 1997;Moldovan et al., 2000;Radulian et al., 2006).These investigations do not consider the influence of the local site responses, but rather average those out.We therefore took a different approach.
As the attenuation relation (Eq.2) does not consider azimuthal differences in attenuation, a new empirical approach is chosen to take the directionality of attenuation and local site effects into account.We introduced a factor (Ardeleanu et al., 2005) into Eq.( 2).This new attenuation model for Vrancea intermediate depth earthquakes is given by the relation: For the determination of the following digitized macroseismic data were used: -Vrancea 10 November 1940 (I epic = IX MSK, moment magnitude M w =7.7, h=150 km, Fig. 2): The macroseismic map of Demetrescu (1941) with the intensity distribution in Romania and Bulgaria was digitized and georeferenced by the Institute for Photogrammetry and Remote Sensing (IPRS) of Karlsruhe University.The data of the Republic of Moldova and Ukraine were taken from the map of Sagalova (for details see Bonjer et al., 2008).According to Radu et al. (1990), Radu and Sandi (1992) as well as Atanasiu and Kräutner (1941) reduced intensity values of the original data of Demetrescu and Sagalova were used to meet the MSK-scale and thus enabling comparison with the 1977 and 1986 Vrancea earthquake intensities.
-Vrancea 4 March 1977 (I epic =VIII MSK, M w =7.5, h=94 km, Fig. 3), and 30 August 1986 (I epic =VIII MSK, M w =7.2, h=130 km).Because no listings of the intensities of the Romanian territory were available, the macroseismic maps of Radu and Utale (1989) and Radu et al. (1987) with a scale of 1:1 000 000 were digitized and geo-referenced as well by the IPRS of Karlsruhe University.For the 1977 intensities of the Bulgarian territory an upgraded version (Glavcheva, 1983) of the Grigerova et al. (1978) intensities was used.R. Glavcheva (personal communication) also provided the listing of the Bulgarian intensity data for the Vrancea 1986 event.
In Ardeleanu et al. (2005) the macroseismic data of the 30 May 1990 event were used instead of those of the 1940 event, as we do here.The reason for this exchange is that in the 1990 event two shocks with comparable strength occurred within 14 h.So, the macroseismic effects of both events may be superimposed and therefore inseparable.was calculated for each observation point fixing α=0.001 km −1 .Strong local variations of were avoided by calculating mean values inside grid cells of 0.5 degree in longitude and 0.25 degree in latitude, separately for each event.From these mean grid values the median is taken.Intensities for rectangles without observations were 2-D-interpolated, respectively extrapolated.Figure 4 shows the resulting field.The contribution of the Vrancea intermediate depth source zone to each grid point was computed with the corresponding representative of this point.
Using the assigned value for each point of computation, the seismic hazard of the Vrancea intermediate depth zone (Vi) is calculated in the same way as for the crustal zones.A seismogenic depth of 120 km is assumed for the Vi source zone.Figure 6 shows as an example the seismic hazard from Vrancea intermediate depth earthquakes for the considered region for a recurrence period of 475 years.

Results
The final seismic hazard map for a recurrence period of 475 years in Fig. 7 is compiled by adding the expected value of the annual frequency of a certain intensity for normal depth source zones (Fig. 5) and for the Vrancea intermediate depth zone (Fig. 6) for every grid point before the common probability of exceedance of the final map is computed.The influence of the Vrancea intermediate depth zone on the seismic hazard is predominant for Romania, for Bulgaria it is important only for its northernmost part.
The separate computation and the separate compilation of the maps of the seismic hazard for the crustal and for the intermediate Vrancea earthquakes are not only an interesting feature for discriminating the sources and their contribution to the seismic hazard of a particular site.The separate computation is especially necessary for the creation of site specific response spectra in order to account for the different spectral content: the ground acceleration of the intermediate earthquakes reaches the maximum values in the frequency band of 0.5-3 Hz (Lungu et al., 1999), whereas the crustal earthquakes have their acceleration maxima in a higher frequency band.
The building code EC 8 recommends a recurrence period of 475 years (probability of exceedance of 10% in 50 years) for the design earthquake (Fig. 7).In order to limit the damage of buildings and financial loss in case of weaker earthquakes with higher frequency of occurrence, a second hazard level in EC8 is recommended, corresponding to a recurrence period of 95 years (probability of exceedance of 10% in 10 years).The hazard map for both territories -Bulgaria and Romania -is presented in Fig. 8. Special structures like large dams or high risk facilities are beyond the scope of EC 8, as they require higher safety standards.This would imply longer recurrence periods.In order to give decision makers a first information for regional planning we present in Figs. 9 and 10 the seismic hazard for a recurrence period of 2475 years (2% in 50 years) and 4975 years (1% in 50 years), respectively.These maps of course can not replace site specific expertises.
In comparison to former calculations in Ardeleanu et al. (2007) here the new hazard maps show reduced values in many parts.This is due to the following reasons: the statistics with classes in steps of half intensities (in former calculations half intensities have been counted to the next higher full intensity), enlarging the field to the Bulgarian territory because of new available digitized macroseismic data from Vrancea intermediate depth earthquakes (in former calculations the values for Bulgaria were assumed conservatively), and the downsizing of the Vrancea intermediate depth seismic source zone based on precise hypocenter localization.This in all results in more realistic seismic hazard maps.

Fig. 3 .
Fig. 3. Macroseismic map (in MSK intensity) of the 1977 Vrancea earthquake (4 March 1977) with epicentral intensity I epic =VIII MSK, moment magnitude M w =7.5 and focal depth h=94 km.Intensity values are from Radu and Utale (1989) for the Romanian territory.The Bulgaria data are a refined version of Grigerova et al.(1978).The yellow star represents the rupture initiation (epicenter), the white star the rupture termination (stopping) according toMüller et al. (1978).
3 for the destroying earthquakes of 1940

Fig. 4 .
Fig. 4. The distribution of .The pattern mimics the trends of the general strong decay of the seismic intensity pattern of the big intermediate depth Vrancea earthquakes in NW-SE direction.The green typed values are based on macroseismic observations.Because no detailed digital macroseismic data outside Bulgaria and Romania were available, the values there were fixed by us (red numbers).Therefore the seismic hazard values outside Bulgaria and Romania in Fig. 6-10 may be erroneous.

Fig. 5 .Fig. 6 .Fig. 7 .
Fig. 5. Seismic hazard from source zones of normal depth for a recurrence period of 475 years; colours represent the intensities in MSK-1964 scale.

Fig. 8 .
Fig. 8. Seismic hazard from all source zones (zones of crustal events and Vrancea events of intermediate depth) for a recurrence period of 95 years; colours represent intensities in MSK-1964 scale.The hazard values outside Romania and Bulgaria may be erroneous (see remarks to Fig. 4).

Fig. 9 .
Fig. 9. Seismic hazard from all source zones (zones of crustal events and Vrancea events of intermediate depth) for a recurrence period of 2475 years; colours represent intensities in MSK-1964 scale.The hazard values outside Romania and Bulgaria may be erroneous (see remarks to Fig. 4).

Fig. 10 .
Fig. 10.Seismic hazard from all source zones (zones of crustal events and Vrancea events of intermediate depth) for a recurrence period of 4975 years; colours represent intensities in MSK-1964 scale.The hazard values outside Romania and Bulgaria may be erroneous (see remarks to Fig. 4).

Table 1 .
Parameters of intensity-frequency relations and input parameters for seismic hazard calculation. .