Probabilistic Seismic Hazard Analysis

The reason for writing this article was the information spread by a number of media outlets about the forecast of the “most powerful” earthquake that may occur in the next 30 years in Japan and the Kuril Islands with a probability of up to 40%. Journalists referred to Japanese scientists. We managed to find the original article, from where, apparently, and taken information. It was published on December 20 in Japan News and is now in the paid archive, but we have a wonderful wayback machine resource.

Below we will try to make out what the source is really about, plunging into the wilds of mathematics and the basics of probabilistic seismic hazard analysis. If in a nutshell, Japanese seismologists did not give a forecast of the mega-earthquake in the Kuril Islands, but described the model characteristics of seismic sources, which are taken into account when mapping seismic zoning with the service life of buildings and structures for the next 30 years. We will try to make calculations. Without supercomputers ...

So, the main task of seismic zoning is to determine the most likely seismic effects. The term "most" we decipher below. To begin with, we will write down simple mathematics and introduce some terms that are understandable to most readers, even those who are more or less familiar with probability theory and mathematical statistics. For definiteness, seismic effects will be measured in physical parameters of acceleration. For example, the peak ground acceleration of 0.7g (where g=9.81 /^2 is the acceleration due to gravity) corresponds to the intensity of shocks of about IX points, and IX points are a disaster.

Each earthquake has its own magnitude M. It is a quantity that characterizes the force (gap size, energy, etc.) in the source. Not to be confused with the score! A score is already a manifestation of shakiness on the surface, that is, sensibility. The closer to the epicenter, the more intensively the Earth's surface oscillates, the higher the level of seismic effects, and moving away from the epicenter, the seismic signals attenuate. This is a fundamental property inherent in all signals in a continuous medium. Let's call it the law of attenuation of seismic effects. Obviously, the building is not destroyed by a magnitude, but by peak soil accelerations, which exceed a certain limit set during design and construction (Fig. 1).
')

Figure 1. The aftermath of the 2007 Nevelsk earthquake

The model of peak ground acceleration attenuation (peak ground acceleration) is defined by the function g(m,r) , which determines the dependence of the (natural) logarithm of peak ground acceleration ln⁡ PGA on the event with a magnitude m and at distance r . This function is represented by a regression ratio (Fig. 2), built on the basis of a regional base of accelerograms obtained using a network of accelerographs. Usually it looks like:

$$

$$g (m, r) = ln \ overline {PGA} (m, r) = c_1 + c_2m + c_3ln⁡ (r + c_4) -c_5 r + c_6 F + c_7 S + σ, (1)$$

where c_1, c_2, c_3, c_4, c_5, c_6, c_7 are the regression coefficients, and F and S describe the dependence on the fault type (Fault) and soil characteristics (Soil), respectively.


Figure 2. Attenuation curve of peak ground acceleration with distance for earthquake with magnitude M = 5.8: 1 - measured values, 2 - regional regression, 3 - confidence interval ± σ, 4 - density function of the standard normal distribution of the natural logarithm of peak acceleration at a distance of 10 km from the source, 5 - the level of peak accelerations of 0.7g.

The variability (from event to event and from point to point) observed in the data on strong ground movements is described by the normal distribution of ln ⁡PGA(m,r) at each point (m,r) by means of zero mean and standard error σ (Fig. 2). This should still be treated as an admission, but quite working [Cornell, 1968].

Then, from the assumptions made above regarding the lognormal distribution of the peak acceleration of the soil, it follows that of all the events of magnitude m occurring at a distance r from the point under consideration, the proportion of those events that will not cause seismic effects exceeding the acceleration a is

$$

$$F (\ frac {lna-g (m, r)} \ sigma), (2)$$

where is the distribution function of the standard normal value. The annex is a table with the numerical values ​​of the normal distribution. We borrowed it from [Baker, 2008].

Thus, the probability of not exceeding a given level of seismic effects a from a seismic event with a magnitude m at a distance r will be fully specified by the distribution function of the standard magnitude. We denote this probability as Pr(A≤a│m,r) . Then, in accordance with definition (2), we can write:

$$

$$Pr (A≤a│m, r) = \ int \ limits _ {- \ infty} ^ a \ frac {1} {\ sigma \ sqrt {2 \ pi}} \ exp (- \ frac {1} {2 } (\ frac {lnx-g (m, r)} {\ sigma}) ^ 2) dx. (3)$$

The graphical interpretation of formulas (1) - (3) is shown in fig. 2. It can be seen that the regression describes the measured values ​​of peak soil accelerations. All points refer to one earthquake. More precisely, the regression itself is based on a large set of empirical data of different earthquakes, but for simplicity, graphs are given for earthquakes with a magnitude of M=5.8 . For example, near Nevelsk (Fig. 1) a series of earthquakes of a magnitude close to M~6 occurred, causing destruction in the city.

In fig. σ=0.783 2 shows the scatter of the values ​​of peak acceleration, which determines the standard error σ=0.783 from (1). The natural logarithm of the peak ground acceleration for the points m=5.8 and r=10 is g(m,r)=4.868 (PGA=125 /^2) . Seismic effects with a peak acceleration of 0.7g (687 /^2) correspond to ln⁡ a=6.532 . Then the probability that the seismic effects of A caused by an earthquake with a magnitude of m=5.8 at a distance of r=10 will not exceed 0.7g (IX points) is equal (see Appendix):

$$

$$Pr (A≤0.7g│m = 5.8, r = 10) = (\ frac {6.532-4.868} {0.783}) = 0.983$$

Buildings are not built for endless use. Everything has its own lifetime, and it is determined by the time T


Figure 3. Schematic illustration of seismic sources.

In what case of the impact of an earthquake with a magnitude M_i , the source of which is located at the point k (Fig. 3), does not exceed the magnitude of the acceleration a at a given point in the next period of time T ? The first answer to such a long question is obvious: if there is no earthquake at all in the considered time interval. The second case is less obvious: if one earthquake occurs, but according to the damping law defined in (1), it will not cause the expected seismic effects. In general, during time T , at least two earthquakes or three or even Ns may occur. Let us rephrase what has been said in terms of probability.

We adopt a simplified model in which a certain point generates earthquakes of only one magnitude. Let's look at the table. 1. It shows various outcomes in which seismic effects will not be exceeded: no seismic events (0), one (1), two (1 + 1), three (1 + 1 + 1), etc. The probability of occurrence of one event during the time T is equal to p1 , two - p2 , three - p3 , etc., none - respectively, p0 . Naturally, the considered set of outcomes is complete, i.e. there are no other options:

$$

$$p0 + p1 + p2 + p3 + \ dots = 1. (4)$$

The probability that the oscillations from an earthquake of a given magnitude does not exceed the specified seismic effect at a given point is denoted as Pr , which is completely given by (3). Then the probability that oscillations from two earthquakes do not exceed the specified seismic effect is equal to Pr*Pr , three to Pr*Pr*Pr , etc. In the absence of earthquakes, it is obvious that the probability of not exceeding is equal to 1 . Thus, the total probability P that during the time T seismic effects at a given point will not be exceeded:

$$

$$P = p0 \ ast (1) + p1 \ ast (Pr) + p2 \ ast (Pr \ ast Pr) + p3 \ ast (Pr \ ast Pr \ ast Pr) + \ dots. (5)$$

Table 1. Probabilistic models of the source of earthquakes of a given magnitude at a given point.

Outcomes (number of events)	Probability of occurrence	Probability of not exceeding
0	p0	one
one	p1	Pr
1 + 1	p2	Pr * Pr
1 + 1 + 1	p3	Pr * Pr * Pr
...	...	...

After the expression (5) is completely digested, we can proceed to the consideration of the general case. To do this, we need to set the probability that during the subsequent time T s events will occur, s=0,1,2,…,Ns . Denote the probability of occurrence of s events with a magnitude M_i at some point k in the subsequent time T as P_k(s,M_i,T) .

Then the probability of not exceeding a given level of seismic effects a from a seismic event with a magnitude M_i at point k (Fig. 3) in the subsequent time interval T will be set by analogy with (5):

$$

$$P (A \ leq a \ mid M_i, T, k) = \ sum_ {s = 0} ^ {N_s} P_k (s, M_i, T) [Pr (A \ leq a \ mid M_i, R_k)] ^ S, (6)$$

where R_k is the distance from the earthquake source at point k to the point at which we expect seismic effects.

Further - easier. It is necessary to consider independent implementations of magnitudes M_i . Of course, there are exceptions in high seismology. For example, when a strong earthquake causes a series of aftershocks (or in other words, repeated events with a magnitude slightly less than that of the main event). But for this purpose, the so-called declustering of the earthquake catalog is produced, i.e. removal of aftershocks and other dependent events.

Taking into account the condition of independence of magnitudes from (6), we obtain the probability that a given level of seismic effects a will not exceed a from a series of earthquakes at point k in the subsequent time interval T :

$$$$ display $$ P (A≤a│T, k) = \ prod_ {i = 1} ^ {N_m} Pr (A≤a│M_i, T, k). (7) $$ display $$

It is reasonable to assume that seismic sources are independent, and each of them lives its own life. Again, in high seismology, there are cases when the tectonic activity of one fault causes seismicity to become active on another. Nevertheless, from the assumption of independence of the set N seismic sources, we obtain the probability of not exceeding a given level of seismic effects a in our “waiting” section in the subsequent time interval T :

$$$$ display $$ P (A≤a│T) = \ prod_ {i = 1} ^ {N} \ prod_ {i = 1} ^ {N_m} Pr (A≤a│T, k). (8) $$ display $$

Let us pass from the probability of non-exceedance in (8) to the probability of exceeding:

$$

$$P (A> a│T) = 1-P (A≤a│T). (9)$$

Expression (9) is the base for the production of seismic zoning maps in probabilistic seismic hazard analysis.

Let us briefly P_k (s,M_i ) probability of an earthquake P_k (s,M_i ) . In seismic zoning, an exponential Poisson model is often used, that is, the probability that in the next T years s earthquakes occur is equal to:

$$

$$P_k (s, M_i, T) = \ frac {[λ_k (M_i) T] ^ s exp [-λ_k (M_i) T]} {s!}, (10)$$

where λ_k (M_i ) is the frequency of occurrence of earthquakes with a magnitude M_i at the point k . As it should be, the sum (10) over all s from zero to infinity is one!

Japanese scientists according to the Paleotsunami made the reconstruction of historical earthquakes. As you know, tsunamis are caused by quite strong seismic events, and their splash persists in the geological and morphological history of the coast. According to the estimates of the Japanese in the area of ​​the southern Kuril Islands, the magnitude of the maximum possible earthquake can be at least M=8.8 . As they believe, the events in question do not occur by chance, but periodically. The interseismic interval lies in the range of 340-380 years, i.e. The average recurrence period is 360 years. Then the frequency of occurrence of earthquakes λ will be estimated as 1/360 ^-1 .

Let us again consider a simplified model of a seismic source — at a given point of a given magnitude (Fig. 3). The standard time interval is T=30 . Then it is obvious that λT << 1 . This means that in formula (10) we can restrict ourselves to the “first” probabilities p0 and p1 :

$$$$ display $$ p0 = exp⁡ (-λT), p1 = λTexp⁡ (-λT). (11) $$ display $$

From (11) it follows that at each instant of time the occurrence of an earthquake is equally probable. For example, according to this model, the sensational “powerful” earthquake will occur in the next 30 years with a probability of 0.083 .

The probability of not exceeding the specified seismic effects by the formula (5) will be evaluated as

$$$$ display $$ P = p0 + p1 * Pr = exp⁡ (-λT) + Pr⁡ * λT * exp⁡ (-λT), (12) $$ display $$

The required probability of exceedance according to (9) and (12) is given by

$$$$ display $$ Q = 1-P = 1-exp⁡ (-λT) -Pr⁡ * λT * exp⁡ (-λT). (13) $$ display $$

Expanding (13) in powers of λT and leaving only the first terms of the expansion, we get:

$$$$ display $$ Q (A> a│T) = λT * (1-Pr (A≤a│m, r)). (14) $$ display $$

Thus, under certain assumptions, we obtained a simple formula for estimating the normative seismic effects by the probabilistic method. Remained last detail. This is the law of attenuation of strong ground motion for the Kuril Islands.

For the Kuril Islands, the neighboring region of Japan is similar in geological and tectonic conditions. Here, as in California, the attenuation equations for peak ground accelerations have been developed at the current level. The most unified model that takes into account the division of seismicity into subduction and crust types is the basic model of Shi and Midorikawa [Si, Midorikawa, 1999]. Later modifications of this model are largely related to the refinement of the coefficients after a detailed analysis of the seismic effects of the 2011 Tohoku mega-earthquake (Mw 9.0) .

$$$$ display $$ \ begin {cases} log⁡A = b-log⁡ (X + c) -kX, D≤30 km \\ log⁡A = b + 0.6 log⁡ (1.7D + c) -1.6 log ⁡ (X + c) -kX, D> 30 km, (15) \ end {cases} $$ display $$

Where

$$

$$b = αMw + hD + \ sum d_i S_i + e, (16)$$

A - peak acceleration in /^2 ; D is the “average” depth of the rupture plane (centroid depth) in km; X is the shortest distance to the gap in ; Mw - moment magnitude;


For crustal earthquakes, S_1=1, S_2=0, S_3=0; for interplate earthquakes S_1=0, S_2=1, S_3=0; for intraplate S_1=0, S_2=0, S_3=1.


Figure 4. Model of the seismic source in the southern Kuriles (source)

"Our" mega-earthquake is expected at the junction of two plates, i.e. type of earthquake - interplit. The 2011 Tohoku mega-earthquake, known for the tragedy at the Fukushima nuclear power plant, occurred at a depth of about 11 km. Therefore, in our case, we assume the source depth D=11 , and the magnitude Mw=8.8 (as suggested by the Japanese). In fig. Figure 4 shows the spatial model of the mega-earthquake expected on the Kuriles. So, for example for oh. Shikotan is the shortest distance from the fault plunging in the direction of the Kuril Islands to the center of the island is about 40 km, so we take X=40 . Now that all the input parameters are defined, we define the probability Pr(A≤a│m=8.8,r=40) from (14) according to the damping law (15). To do this, you need to tinker a bit, to bring the expression (15) to the natural logarithm, etc. As a result, we obtain: g(m,r)=7.22 , σ_ln⁡ PGA =0.62 . Below, in the table. 2, the calculated probability values ​​are presented.



The purpose of this analysis is to estimate the probability of exceeding the ground movement level, and the main result is to determine the dependence of the excess probability on the level of movement, which is called the hazard curve. The hazard curve is clearly shown in Fig. 5 according to the data from table. 2. From fig. 5, we see that the probability of exceeding 0.06 over the next 30 years corresponds to a seismic impact level of about 0.9g. This is more than 9 points. Such an estimate for the Kuril Islands is a natural and expected wave, if we recall that here in 1994 there was a major earthquake, called the Shikotan . It was accompanied by destruction of buildings, death of people, tsunami waves and numerous landslides. The intensity of tremors on about. Shikotan ranged from VIII to IX on the MSK-64 scale.




Figure 5. Hazard Curves for Fr. Shikotan on the model of a seismic source in the southern Kuriles according to this source .

In addition to the Poisson model, there are models of the occurrence of earthquakes with “memory”, which take into account the history of previous seismicity. These include the Brownian model (Brownian Passage Time). It is used by Japanese colleagues when building seismic zoning maps for well-studied territories. According to their estimates, the probability of mega-earthquakes in the next 30 years in the southern Kurils was from 0.07 to 0.4 (Fig. 4).

Imagine that a maximum earthquake occurs strictly periodically - once every 1000 years. And here it happened, say, 20 years ago, and we are going to build residential buildings in this place, which will serve the next 50 years. During this time, strong events will not happen anymore, as they were just now, and the next one will be repeated no earlier than in 980 years. If this fact of periodicity is established, then the Brown model is used. Then we will get seismic zoning maps with lower values ​​of the calculated intensity of seismic impacts, and this will reduce the cost of construction at the output. What model to choose, decides not one specialist, but a group of experts.

To calculate the hazard curve according to the Brownian model, more time and materials are required. However, we will accept for simplicity that the probability of an earthquake occurring in the next 30 years is known to us and it is equal to p1, according to a report from Japanese colleagues. Then the probability that no events will happen in the next 30 years will take p0 = 1-p1. Finally, the probability of exceeding will be set as

$$$$ display $$ Q (A> a│T) = 1-P = 1-p0-Pr⁡ * p1 = p1 * (1-Pr (A≤a│m, r)). (17) $$ display $$

For the maximum probability of a mega-earthquake of 0.4, we constructed a hazard curve (Fig. 5). It shows the exorbitant values ​​of the design seismic effects (Fig. 5), more than X points. In reality, we did not take into account many points concerning ground conditions, the nature of the attenuation of seismic effects from mega-earthquakes, etc.It is known that seismic effects from earthquakes with a magnitude are Mw=9.0not much more than those caused by events with Mw=8.3. This means that we clearly overestimated the seismic effects by the formula (15). Nevertheless, the presented calculation allows you to feel the scale of the problem.

Last thing.Why did Japanese seismologists give a probabilistic model of seismic sources for the next 30 years? In Japanese standards, the operational terms of the “life” of buildings and structures are 30 and 50 years, in Russian construction rules - 50 years. The technology of seismic zoning around the world (including in our country) comes down to an assessment of seismic effects, which will be exceeded for 30 or 50 years with a given probability. For objects of normal responsibility (and these are our houses and offices with you), the building codes set the probability of exceeding 0.1 for 50 years. For 30 years, the probability is 0.06 (this is Japan and some other countries).

Thus, Japanese seismologists have adopted a conservative estimate of the seismic model, which will be used when updating seismic zoning maps. This underlines the high responsibility and culture of seismological studies of Japanese scientists. It also suggests that the Earth knows no boundaries, and the study of natural phenomena requires the joint efforts of scientists from different countries.

The analysis was performed by Alexey Konovalov (a.konovalov@geophystech.ru), deputy director for research at GEOFIZTECH LLC.