A method for Bayesian estimation of the probability of local intensity for some cities in Japan

Seismic hazard in terms of probability of exceedance of a given intensity in a given time span, was assessed for 12 sites in Japan. The method does not use any attenuation law. Instead, the dependence of local intensity on epicentral intensity I0 is calculated directly from the data, using a Bayesian model. According to this model (Meroni et al., 1994), local intensity follows the binomial distribution with parameters (I0, p). The parameter p is considered as a random variable following the Beta distribution. This manner of Bayesian estimates of p are assessed for various values of epicentral intensity and epicentral distance. In order to apply this model for the assessment of seismic hazard, the area under consideration is divided into seismic sources (zones) of known seismicity. The contribution of each source on the seismic hazard at every site is calculated according to the Bayesian model and the result is the combined effect of all the sources. High probabilities of exceedance were calculated for the sites that are in the central part of the country, with hazard decreasing slightly towards the north and the south parts.


Introduction
The distribution of seismic intensity is generally influenced by major geological and tectonic features and, on a local scale, by local geological conditions.Factors such as fault direction, depth of focus, local geological and topographical conditions, small scale irregularities, are among others of primary importance.

Key words Bayesian estimation -local intensityrandom variable p -probabilities of exceedance -Japan
Many relationships have been proposed by different authors to model the variations of macroseismic intensity versus distance (Algermissen et al., 1976;Anderson, 1978;Chandra, 1979;Ambraseys, 1985;Grandori et al., 1991;Papazachos et al., 1993;Musson and Winter, 1997;Slejko et al., 1998) among others.Their leit-motiv is that macroseismic data could, in some way, be used like instrumental ones for determining earthquake source parameters and properties of wave propagation Utsu (1988) derived a relation between the seismic intensity near the epicenter, the focal depth and the magnitude of earthquakes which is of the form (1.1) this is an empirical formula and holds for seismic intensity (which is in JMA-Japanese-scale) I 0 > 0, for M ranged between 2.0 and 8.0 and focal depth h from 3 up to 100 km.Utsu noted that the magnitude used was the M (of JMA) which is roughly equal to M S in the interval 5.0 to 8.0.The eq. (1.1) was determined using 1114 data (earthquakes).The decay of intensity with distance and its spatial distribution are performed when the observed intensities of an earthquake are a sufficient number of points.A number of attenuation laws were constructed for numerous regions of the world, and their parameters have been evaluated using simple radiation model (s) involving a point source and by applying a least squares technique to the radii D i of the isoseismals.In order to model the uncertainty in the data (Meroni et al., 1991) the propagation itself of the intensity, from the epicenter to a site, has been considered as a process which is affected by randomness due to the qualitative character of the seismic intensity, as well as to the complex relationships between intensity decay and local seismotectonic features.In order to make this aspect clear the normalized intensity decay D i is considered as a random variable following the Beta distribution given the epicentral intensity I 0 and the distance R from the epicenter to the site of interest.
The paper confines itself to the estimation of the probability distribution of the intensity I S at some cities of Japan, only considering the observed intensity points (fig. 1) of a sufficient number of earthquakes (fig.2) belonging to a zone which has the same attenuation.No complex attenuation laws are needed because in the study the estimates are based only on the data.Table I shows the events used with the available observed intensity data points occurring in different cities which are within the zones of Japan.The number of the data points for each intensity level is listed in the right part of the table.

Method applied
The observed intensity at the site, I S , is considered as a random variable having the binomial distribution.In our case, the binomial distribution simulates values of probability (in the same way that an exponential or power function simulates the values of intensity in an attenuation law, when using a non-probabilistic approach).That is, if I 0 I 0 and epicentral distance r.It is sensible to discretize the values of epicentral distance, r j and round all epicentral distances r to the closest r j .The values of epicentral distance increase with a step of 10 km.
In order to apply a Bayesian model the parameter p must be considered as a random variable, as well.A prior Beta distribution is required for p and is expressed by where α and β are the Beta distribution parameters, Γ is the Gamma function and x is the integration variable.Meroni et al. (1994) also proposed that the prior expected average of α and β for epicentral intensity I 0 = i 0 and epicentral distance r = r j are (2.3) where and C is a constant such that the prior estimate of p, is almost equal to 0.99.The physical meaning of this, is that for the smallest value of j ( j = 1) which corresponds to the smallest epicentral distance (r 1 ), it is almost certain that the intensity will be equal to the epicentral intensity.
Equations (2.2), (2.3) and (2.4) combined, have as a result that the prior estimator of p is equal to (2.5)It was found that the posterior values of p are not sensible to changes in the prior values of α and β, so we adopted the proposed prior values.It follows that the posterior mean of p is (2.6)where α and β are the parameters of the Beta distribution of p, i S is the n-th intensity point at an epicentral distance r j and N j is the number of observed intensities caused by earthquakes having an epicentral intensity of I 0 and an epicentral distance r j .This is the computed p value, used in the calculation of the probability of exceedance of given intensity.Using eq.(2.6) we calculated values of p for epicentral intensities 3-6 (JMA) and epicentral distances 0-500 km.Following Meroni et al. (1994) if there were no data for a value of r j , we assigned to it the same p with r j-1 .The prior and the posterior estimates of the parameter p, for the epicentral distance 0-500 km, are shown in fig.3a and fig.3b, respectively.For distances greater than 500 km there appears to be no dependence of p on distance (fig.4), mainly because of lack of data.However a large value of intensity in the long distance exists in fig. 4.This can due to an error in the estimation of intensity for any reason (like the wrong determination of shock's epicenter or/ and its corresponding magnitudes).
An inspection to figs.3a,b and 4 illustrates that the posterior estimates of p do not decrease monotonically with increasing distance from the source.A relevant Bayesian theory was given by Rhoades and Evison (1993, eq.(2.10)).The main difference from the present study is that instead of a prior Beta distribution for a single value p, the above mentioned authors assume a prior Dirichlet (multivariate Beta) distribution for a collection of ordered binomial parameters.An alternative approach would be to smooth the posterior estimates of p over radii, for a given value of I 0 .For this purpose we applied smoothing by 7-point moving average.An example prepared for Tokyo where we estimated the probability of exceedance used for this purpose the alternative smoothing procedure (in table III, as Tokyo A).No significant difference from the original results observed.The procedure suggested by Rhoades and Evison (1993), would be very useful if we focused to the attenuation and not to hazard assessment.In our case the parameter p in not strictly decreasing with distance, but its mean value shows such trend.In order to base the results on reliable data only, it was decided to choose 500 km as the upper limit for epicentral distance.In fig. 5 we can see the smoothed posterior value of p for the chosen distance of 500 km.  the computations the a and b values very recently stimated by Koravos (2000) for the zones in which Japan and the surrounding area are divided.In table II the values of a and b parameters used, are listed.
In order to take into account the uncertainties of eq.(1.1), it was assumed that its values have a standard error of σ = 0.5.So, instead of obtaining a single value for M, we calculated a normal distribution with a mean value of M and standard deviation σ.Applying this, the probability P(I 0 ≥ i 0 ) is the one obtained by the average value of eq.(2.8), weighted over all possible values of M using a weighing factor equal to N(M,σ 2 ), where N is the probability function of the normal distribution, M is the magnitude given by (eq.(1.1)) and σ is the standard error (2.9) where P(I 0 ≥ i 0 | m) is the probability of exceedance of epicentral intensity I 0 given that the magnitude is m and N(m; M, σ) is the probability density function of the Normal distribution at m.
In comparison with the method proposed by Meroni et al. (1994), we showed an alternative approach.These authors utilised the probability distribution of the epicentral intensity.Contrary to the aforementioned method we express the probability distribution of I 0 through the probability distribution of M i (eq.(2.8)) in combination with Utsu's magnitude-epicentral intensity relation (eq.(1.1)).In addition to this we assumed that there is an uncertainty in this relation, which is incorporated in our method with the use of eq.(2.9).
The second factor of eq.(2.7), is calculated by eq.(2.1) with the proper value of p for the specific epicentral distance and epicentral intensity.
For every subdivision of every source, the probability in question is the sum of the probabilities that correspond to epicentral intensities of I 0 = I S , I 0 = I S + 1, I 0 = I S + 2 and so on, up to the maximum value of I 0 , which is 6 for the Japanese intensity scale (2.10) In order to assess the seismic hazard in a specific site, in terms of probability of exceedance of a given intensity in a given time span, the following steps must be taken.First, for every subdivision of every source and for every possible value of I 0 , the probability of exceedance of I S in the given site is evaluated (2.7) where P(I 0 = i 0 ) is the probability of occurrence of an earthquake with an epicentral intensity of I 0 , in the given subdivision of the given zone within the given time span and is the probability of an earthquake with an epicentral intensity of I 0 , to cause an intensity of I S or greater at the site.
The first factor of eq.(2.7) is calculated indirectly.The magnitude of an earthquake that is expected to have an epicentral intensity of I 0 , is given by eq.(1.1).A mean value of h equal to 15 km was calculated from the data.Then using the Gutenberg-Richter law parameters a and b (1944) for the specific source, the probability of exceedance of the magnitude M is calculated (2.8) where M 0 = M (i 0 ) see eq. (1.1) and then it is normalized for an area equal to the area of the subdivision.In the present study we adopted for  Table III.Probabilities of exceedance for the site intensities 6.0, 5.0, 4.0 and 3.0 in 1, 2, 5, 10, 20, 25, 50, 100, 200 and 500 years for the 12 Japanese cities examined.The epicentral distances vary from 0 to 500 km.An example prepared for Tokyo (see Tokyo A) used for this purpose the smoothing by 7-point moving average.
The contribution, in terms of probability of a whole source is the sum of the probabilities of all the subdivisions of the source.Finally, the contributions of all the sources are combined using the relation (2.11) where P j is taken as a computation from eq. (2.10).
The seismicity of the examined area is of the highest in the world and observed to occur in the front and the back-arc area where great interplate and intaplate earthquakes are generated, respectively.In both areas, high risk seismic sources were created because of their proximity to major centers of dense population and heavy industry.The probability of exceedance of a given intensity in a given time span was assessed for 12 sites in Japan (fig.6).These sites belong both to the front and the back-arc parts of Japan and they are found within zones of the area defined by Papazachos et al. (1994) and Matsuda (1990), respectively.These sites are the following cities (from north to south): Sapporo, Kushiro, Hachinohe, Sendai, Niigata, Hitachi, Tokyo, Gifu, Osaka, Hiroshima, Kochi and Fukuoka.

Results and discussion
Figure 6 shows the seismogenic sources of Japan (Matsuda, 1990;Papazachos et al., 1994) and the 12 examined sites (cities).The method effectively estimates the probability of the local intensity I S through the Bayes statistics.There is no need for an attenuation law, but the seismic hazard can be evaluated directly from the observed data points.The obtained results can be considered as a measure (under Poissonian hypothesis) of the uncertainty in seismic hazard evaluation.Looking at fig. 7, presenting the probability of exceedance of various intensities I S (VI, V, IV and III) for the next 50, 100, 200 and 500 years, we can observed that Kushiro is most likely to experience an intensity VI degrees of the Japanese scale.Hitachi and Tokyo follow, while Fukuoka has the lowest probability of exceedance of an intensity VI in the next 50 years.The seismic hazard of Japan is illustrated in fig.8, where the probabilities of exceedance in the next 50 years (475-years mean return period) of an intensity VI for the twelve examined sites-cities are presented by black circles.This plotted curve taken into account the cities from north towards south and thus fig.8 is an interesting picture of the fluctuations of the seismic hazard profile of Japan.We can observe that the general trend of the hazard profile decreases southwards.Another interesting observation is that almost all the sites-cities in the back arc area are characterized by low probability of exceedance.These cities are: Sapporo, Niigata, Gifu, Hiroshima and Fukuoka.But there are also cities which belong to the front area of Japan, where take place the underthrusting process, like Hochinohe and Sendai which are cited in the inland part of Japan away from the centers of the very large earthquakes and thus their values on the hazard curve are low.The highest values appeared in the cities of Kushiro, Fig. 6.The seismic zones into which Japan was divided (Matsuda, 1990) and the 12 sites-cities for which the seismic hazard is estimated.Tokyo and Hitachi.The problem becomes greater for Tokyo because it is either a high dense populated city and moreover there are centers of high risk around (heavy industry, electricity power and chemical plants, etc).
The upper part of fig.8 demonstrates the seismic hazard profile of Japan, where the probabilities of exceedance in the next 200 years (1000years mean return period) of an intensity VI for the twelve examined sites-cities are presented by black triangles.It is expected to follow almost the same pattern as the hazard profile for 50 years, admittedly with the exception of Niigata.Niigata, number 5 (fig.8), shows relatively higher hazard, in comparison with the previous profile it belongs to those sites of low hazard, and this can be due to the different time period examined.Another critical point is that Tokyo still dominates in the range of the high probabilities of exceedance (the seismic hazard profile).
The majority of earthquakes, particularly the large ones (M ≥ 8.0) occur along the Japan trench and Nankai and Sagami troughs.Earthquakes along the plate boundaries generally show lowangle thrust mechanisms and are interpreted to denote that the Pacific and Philippine sea plates are subducting beneath the Eurasian plate (Ando, 1975;Scholz and Kato, 1978).Earthquakes are also occur in the back-arc area.These are referred to as intraplate seismicity (Wesnouski et al., 1982).The largest intraplate events are generally about one magnitude unit less (7 < M < 8) than the great interplate earthquakes, although they are a major danger because of their proximity to centers of high population.
An integration of geological and seismological data for probabilistic seismic hazard in Japan was carried out by Wesnousky et al. (1984).In particular, he estimated the probabilities of the ground shaking of JMA intensity ≥ V during the next 20, 50, 100 and 200 years, induced by either interplate and intraplate earthquakes.For the 50 year period, he suggested Nankai trough and portions of the Japan trench adjacent to east coast on northeast Honshu will experience major events.The 200 years map projects the probable occurrence of a major earthquake along Sagami trough.The cities of Tokyo and Hitachi (point 6 and 7 in fig.8) are very close to Sagami trough and both can be considered sites of high seismic risk.Especially for the back-arc area and for the next 200 years, Wesnousky et al. (1984) estimated high probabilities in the north and central part of Honshu.The seismogenic sources J-8 and J-9 (fig.6) are parts of this high seismic risk area, where the cities of Sapporo and Niigata are located.In the present study, both of these cities show a relatively high probability of exceedance for the next time span of 200 years.The regions of high probability expand as longer time periods are sampled.The results obtained by Wesnousky et al. (1984) in general are not directly comparable with ours, mainly because we search the ground shaking in sites (cities), whereas Wesnousky and his colleagues covered the whole territory of Japan.Despite this, we can observe that few of the cities examined in this study which present high seismic hazard parameter are within the areas of high risk according to Wesnousky et al. (1984).Rikitake (1991) estimated the probability of exceedance for intensities greater than or equal to V and VI of the JMA within the next 10 years in the broad area of Tokyo.His results for intensity ≥ V are 0.40 while ours is 0.23, and for intensity ≥ VI he found 0.049 while our estimation is 0.043.
The plots of fig.7 in comparison with table III, give useful information for the seismic risk assessment in Japan for Earthquake Resistant Design (ERD).

Fig. 1 .
Fig. 1.Sites-points where the observed intensity is reported.

Fig. 2 .
Fig. 2. The examined earthquakes which caused the intensities studied in the present work.

FigFig. 4 .Fig. 5 .
Fig 3a,b.a) The prior estimates of the p parameter given the epicentral distances 0-500 km and I 0 = VI, V, IV and III, and (b) the posterior estimates of the parameter p given the epicentral distances 0-500 km and I 0 = VI, V, IV and III.

Fig. 7 .
Fig. 7.The probability of exceedance of various intensities I 0 = VI, V, IV and III in the next 50, 100, 200 and 500 years for the 12 examined sites-cities of Japan.

Fig. 8 .
Fig. 8. Seismic hazard curves for the next 50 years (475-years mean return period) and for the next 200 years (1000-years mean return period) presented by black circle and black triangles, respectively for the 12 examined sites-cities of Japan.The name of the cities and their corresponding point number on the curves are: 1 -Sapporo; probability of exceedance of intensity I in 500 years probability of exceedance of intensity I in 100 years probability of exceedance of intensity I in 50 years probability of exceedance of intensity I in 200

Table I .
Events used with the available observed intensity data points occurring in different cities within the zones of Japan.The number of data points for each intensity level is listed in the right part of the table.