**2. Construction of generalized vicinity of strong earthquake**

We have used the Harvard worldwide seismic moment catalog for 1976–2005, and the USGS/NEIC catalog for 1968–2007. In both cases only shallow earthquakes with depth H < 70 km were examined. Two subsets of data can be used, first one includes all earthquakes from the catalog and the second includes stronger earthquakes that are only completely reported. Below we present the results from processing of the Harvard catalog using the first subset of data (all reported events) and the results for the USGS/NEIC catalog using only completely reported events. In the latter case the events with magnitude M ≥ 4.7 were used, a total number of events was 97615. A similar cutoff for the Harvard catalog would reduce the available data too much to get statistically robust results.

Both used data sets were searched for events falling into the space–time domains surrounding the source zones of large (M7+) earthquakes, with due account for the seismic moment in the Harvard catalog and the maximum magnitude for the USGS/NEIC catalog. A generalized vicinity of large earthquake is understood as a set of events falling into the zone of influence of any of these strong earthquakes. The zones of influence were defined as following, see also (Rodkin, 2008) for the details. Spatial dimensions of the zones of influence for earthquakes of different magnitudes were calculated from the approximate relationship (Sobolev & Ponomarev, 2003) between typical source size *L* and earthquake magnitude M:

$$L \text{ (km)} = 10^{0.5M - 1.9} \text{.} \tag{1}$$

In the examination below the earthquakes located at distances within 7×*L* from the epicenter of the given strong earthquake were taken into account.

For constructing a time vicinity of strong earthquake we used the conclusion that duration of a failure cycle weakly depends on earthquake magnitude (Smirnov, 2003). Hence the simple epoch superposition method can be used for comparing the time vicinities of earthquakes with close magnitudes. At the figures below all earthquakes located in the area 7×*L* of the corresponding strong event were taken into account. This choice allows the most complete use of available data. Negative consequences of this choice are a lower statistical significance at the edges of the time interval because of shortage of data there, and a false effect of a systematic growth of a number of earthquakes towards the centre of the used time interval. However these errors can be taken into account, so they do not distort the results.

The generalized vicinity of large earthquake which was constructed contained more than 60000 earthquakes for the Harvard catalog and more than 300000 earthquakes for the USGS/NEIC catalog. Such a big number of events resulted from the fact that one and the same earthquake can belong to the space–time vicinities of different strong earthquakes.

Current State of Art in Earthquake Prediction, Typical Precursors and

Experience in Earthquake Forecasting at Sakhalin Island and Surrounding Areas 47

Fig. 1. The fore- (upper panel) and aftershock (lower panel) sequences in the generalized vicinity of strong earthquake, events flow is given in a number of events per day, zero time

Note that as it can be seen in Fig. 1, the mean duration of the foreshock process is significantly shorter than that of the aftershocks, and the rate of increase for foreshocks toward the moment of main shock occurrence is noticeably slower than the decay of the aftershock rate. From this it follows that the typical maximum rate of foreshock sequence is an order smaller than the maximum rate of aftershock process. This result is intimately related to the problem of predictability of large earthquakes; the prediction problem would have been solved already, if the rate of foreshocks would be equal to the rate of aftershocks. It seems important to note that upon closer examination the seismicity increase in the generalized vicinity of large earthquakes is not confined to the foreshock and aftershock cascades. Essentially weaker but quite noticeable increase above the mean rate level occurs

corresponds to the moment of occurrence of the generalized main shock.

Such an increase in a number of events has considerably enhanced the possibility of statistical examination.

Time and space position of each earthquake falling into the generalized vicinity of large earthquake is characterized by the time shift from the origin time of the corresponding strong earthquake and by the distance from the epicenter of this main event (norm to the source size of this main event). Both catalogs (USGS/NEIC and Harvard) were used for examination of the relative space-time density of earthquakes. The Harvard catalog was used in this paper mostly for the verification of results, which were obtained from examination of the USGS/NEIC catalog.

#### **3. Regularities in rate of fore- and aftershock cascades**

The most well known feature of seismic behavior occurring in the vicinities of large earthquakes is the existence of aftershock and foreshock power-law cascades. Figures 1*a* and 1*b* show the foreshock and aftershock sequences in the generalized vicinity of strong earthquake, which were obtained from USGS/NEIC data (similar results were obtained from examination of the Harvard catalog). The earthquakes rate is presented by time density of group of earthquakes consisting of subsequent 50 events taken with step 25 events (rate of events is given in n/day for convenience).

As can be seen in Fig. 1, the evolution of foreshocks and aftershock sequences well correlates with a power law. The Omori law (Utsu et al., 1995; Sobolev, 2003) is known to be a good fit to the aftershock rate:

$$n \sim 1/\ (c+t)\text{-}\mathfrak{v},\tag{2}$$

where *n* is the rate of aftershock occurrences, *t* is the time interval after the main shock occurrence, *c* – parameter fitting the rate of earthquakes in the closest vicinity of the main shock, and *p* is the parameter of the Omori law. The Omori law (2) is a good fit for the interval until one hundred days or somewhat later after the main shock occurrence.

The foreshock cascade occurring before the main shock time can be described in a similar manner; in this case *t* is the time before the main shock origin, and c = 0. The foreshock cascade was found to be quite noticeable in the generalized vicinity 10-20 days before the main shock occurrence (Fig. 1).

Of special interest is the deviation of the aftershock rate from the power law during the first hours after the main shock occurrence. The deficit of earlier aftershocks described by parameter c in (2) is explained sometimes by difficulty in recording all of too numerous aftershocks occurring immediately after a large earthquake. However, this factor is hardly capable of providing a full explanation of the phenomenon (Lennartz et al., 2008; Shebalin, 2006). The deviation from the power law toward lower rates of events during a few first hours of the aftershock sequence can be seen clearly in Fig. 1b; the rate of aftershocks reaches the values obeying the power law only 2–3 hours after the strong earthquake. At that time the mean rate of earthquakes with M ≥ 4.7 occurring in the vicinity of a mean large earthquake (but not in the generalized vicinity of large earthquake) is a little above one event per hour. Such rate can not cause any problem in events recording. Thus, the effect of a lower rate of earlier aftershocks probably has a physical nature. This conclusion is similar with those presented in (Lennartz et al., 2008; Lindman et al., 2010; Shebalin, 2006).

46 Earthquake Research and Analysis – Statistical Studies, Observations and Planning

Such an increase in a number of events has considerably enhanced the possibility of

Time and space position of each earthquake falling into the generalized vicinity of large earthquake is characterized by the time shift from the origin time of the corresponding strong earthquake and by the distance from the epicenter of this main event (norm to the source size of this main event). Both catalogs (USGS/NEIC and Harvard) were used for examination of the relative space-time density of earthquakes. The Harvard catalog was used in this paper mostly for the verification of results, which were obtained from

The most well known feature of seismic behavior occurring in the vicinities of large earthquakes is the existence of aftershock and foreshock power-law cascades. Figures 1*a* and 1*b* show the foreshock and aftershock sequences in the generalized vicinity of strong earthquake, which were obtained from USGS/NEIC data (similar results were obtained from examination of the Harvard catalog). The earthquakes rate is presented by time density of group of earthquakes consisting of subsequent 50 events taken with step 25 events (rate of

As can be seen in Fig. 1, the evolution of foreshocks and aftershock sequences well correlates with a power law. The Omori law (Utsu et al., 1995; Sobolev, 2003) is known to be a good fit

 *n* ~ 1/ (*c* + *t*)–*p*, (2) where *n* is the rate of aftershock occurrences, *t* is the time interval after the main shock occurrence, *c* – parameter fitting the rate of earthquakes in the closest vicinity of the main shock, and *p* is the parameter of the Omori law. The Omori law (2) is a good fit for the

The foreshock cascade occurring before the main shock time can be described in a similar manner; in this case *t* is the time before the main shock origin, and c = 0. The foreshock cascade was found to be quite noticeable in the generalized vicinity 10-20 days before the

Of special interest is the deviation of the aftershock rate from the power law during the first hours after the main shock occurrence. The deficit of earlier aftershocks described by parameter c in (2) is explained sometimes by difficulty in recording all of too numerous aftershocks occurring immediately after a large earthquake. However, this factor is hardly capable of providing a full explanation of the phenomenon (Lennartz et al., 2008; Shebalin, 2006). The deviation from the power law toward lower rates of events during a few first hours of the aftershock sequence can be seen clearly in Fig. 1b; the rate of aftershocks reaches the values obeying the power law only 2–3 hours after the strong earthquake. At that time the mean rate of earthquakes with M ≥ 4.7 occurring in the vicinity of a mean large earthquake (but not in the generalized vicinity of large earthquake) is a little above one event per hour. Such rate can not cause any problem in events recording. Thus, the effect of a lower rate of earlier aftershocks probably has a physical nature. This conclusion is similar

interval until one hundred days or somewhat later after the main shock occurrence.

with those presented in (Lennartz et al., 2008; Lindman et al., 2010; Shebalin, 2006).

statistical examination.

examination of the USGS/NEIC catalog.

events is given in n/day for convenience).

to the aftershock rate:

main shock occurrence (Fig. 1).

**3. Regularities in rate of fore- and aftershock cascades** 

Fig. 1. The fore- (upper panel) and aftershock (lower panel) sequences in the generalized vicinity of strong earthquake, events flow is given in a number of events per day, zero time corresponds to the moment of occurrence of the generalized main shock.

Note that as it can be seen in Fig. 1, the mean duration of the foreshock process is significantly shorter than that of the aftershocks, and the rate of increase for foreshocks toward the moment of main shock occurrence is noticeably slower than the decay of the aftershock rate. From this it follows that the typical maximum rate of foreshock sequence is an order smaller than the maximum rate of aftershock process. This result is intimately related to the problem of predictability of large earthquakes; the prediction problem would have been solved already, if the rate of foreshocks would be equal to the rate of aftershocks. It seems important to note that upon closer examination the seismicity increase in the generalized vicinity of large earthquakes is not confined to the foreshock and aftershock cascades. Essentially weaker but quite noticeable increase above the mean rate level occurs

Current State of Art in Earthquake Prediction, Typical Precursors and

value is calculated from

Experience in Earthquake Forecasting at Sakhalin Island and Surrounding Areas 49

used to examine the change in b-value in the generalized vicinity of strong earthquake (the similar results were obtained from examination of the Harvard catalog). The maximum likelihood method was used for the *b*-values estimation (Utsu, 1965). By this method the *b*-

where *Mav* is the average magnitude for each subset of data and *Mc* is the lower magnitude limit used in the analysis, here *M*c=4.7. Discreteness of magnitude values because of

Fig. 3. The change of mean b-values in the generalized vicinity of strong earthquake. The

values obtained for 50 events groups are given by dots. Panels as in Fig. 1.

*b* = *lg*(*e*) /(*Mav – Mc*) (3)

in the time interval about ±100-300 days around the main shock date. The analysis of the Harvard seismic moment catalog gives similar results; however, these data testify for a broader area of seismicity increase, roughly within ±500 days around the main shock date. This type of long-term pre- and post-shock seismic activity agrees with the suggestion that the final time interval of strong earthquake preparation prolongs a few years.

We now characterize the changes in the seismicity increase as functions of the distance to the main shock epicenter. The distances are compared in units of the magnitude-dependent main shock source dimension *L* from equation (1). Fig. 2 shows the distance–time diagram of rate of a number of events in the vicinity of the main shock. The horizontal axis indicates the time (in days) from the main shock occurrence time; an analogue of longtime scale is used near the main shock occurrence moment. The vertical axis indicates the distance from the main shock epicenter in units of earthquake source size L. Events' rate is given in logarithmic scale, ln(*n*), where *n* is the number of events in a cell of the distance–time diagram.

Fig. 2. Spatial-temporal change of number of earthquakes lg (density of events per day) in the generalized vicinity of strong earthquake, distance Rnorm from the main shock epicenter is given in norm source size units.

As can be seen in Fig. 2, the rate of earthquakes in the vicinity of the main shock begins to increase one-three hundred days before the main shock, and this activity increase accelerates toward the moment of the main shock. The increase in seismic activity occurs at a distance of about three source sizes *L* from the main shock. This estimate of the radius of influence is in agreement with the size of the areas where predictive functions are usually estimated in earthquake prediction algorithms (Kossobokov, 2005; Shebalin, 2006). Seismic activity outside the zone of 3–4 earthquake source sizes decreases in the close time vicinity of strong earthquake. This feature can result from softening in the source of ongoing strong earthquake. Some other indications of strength decrease in the strong earthquake vicinity are presented in (Rodkin, 2008). In this case one can expect the strain rearrangement from the outer "rigid" region into the inner "soft" zone where the strong earthquake is about to occur. As a result of this rearrangement, the probability of earthquake occurrence in the outer "rigid" zone of the future rupture would become somewhat lower.

The decrease in *b*-value is known to be used as an indicator of an increase in probability of a strong earthquake occurrence (Shebalin, 2006; Zavyalov, 2006). Catalog USGC/NEIS was

in the time interval about ±100-300 days around the main shock date. The analysis of the Harvard seismic moment catalog gives similar results; however, these data testify for a broader area of seismicity increase, roughly within ±500 days around the main shock date. This type of long-term pre- and post-shock seismic activity agrees with the suggestion that

We now characterize the changes in the seismicity increase as functions of the distance to the main shock epicenter. The distances are compared in units of the magnitude-dependent main shock source dimension *L* from equation (1). Fig. 2 shows the distance–time diagram of rate of a number of events in the vicinity of the main shock. The horizontal axis indicates the time (in days) from the main shock occurrence time; an analogue of longtime scale is used near the main shock occurrence moment. The vertical axis indicates the distance from the main shock epicenter in units of earthquake source size L. Events' rate is given in logarithmic scale, ln(*n*), where *n* is the number of events in a cell of the distance–time

Fig. 2. Spatial-temporal change of number of earthquakes lg (density of events per day) in the generalized vicinity of strong earthquake, distance Rnorm from the main shock epicenter

As can be seen in Fig. 2, the rate of earthquakes in the vicinity of the main shock begins to increase one-three hundred days before the main shock, and this activity increase accelerates toward the moment of the main shock. The increase in seismic activity occurs at a distance of about three source sizes *L* from the main shock. This estimate of the radius of influence is in agreement with the size of the areas where predictive functions are usually estimated in earthquake prediction algorithms (Kossobokov, 2005; Shebalin, 2006). Seismic activity outside the zone of 3–4 earthquake source sizes decreases in the close time vicinity of strong earthquake. This feature can result from softening in the source of ongoing strong earthquake. Some other indications of strength decrease in the strong earthquake vicinity are presented in (Rodkin, 2008). In this case one can expect the strain rearrangement from the outer "rigid" region into the inner "soft" zone where the strong earthquake is about to occur. As a result of this rearrangement, the probability of earthquake occurrence in the

The decrease in *b*-value is known to be used as an indicator of an increase in probability of a strong earthquake occurrence (Shebalin, 2006; Zavyalov, 2006). Catalog USGC/NEIS was

outer "rigid" zone of the future rupture would become somewhat lower.

the final time interval of strong earthquake preparation prolongs a few years.

diagram.

is given in norm source size units.

used to examine the change in b-value in the generalized vicinity of strong earthquake (the similar results were obtained from examination of the Harvard catalog). The maximum likelihood method was used for the *b*-values estimation (Utsu, 1965). By this method the *b*value is calculated from

$$b = \lg(e) \;/\; (M\_{av} - M\_c) \tag{3}$$

where *Mav* is the average magnitude for each subset of data and *Mc* is the lower magnitude limit used in the analysis, here *M*c=4.7. Discreteness of magnitude values because of

Fig. 3. The change of mean b-values in the generalized vicinity of strong earthquake. The values obtained for 50 events groups are given by dots. Panels as in Fig. 1.

Current State of Art in Earthquake Prediction, Typical Precursors and

earthquakes …, 1976).

forecasting" is used.

**4.1.2 Methodology** 

**earthquake (the South Kuril region) 4.1.1 Seismic region and data** 

inefficient.

Experience in Earthquake Forecasting at Sakhalin Island and Surrounding Areas 51

scientists should also assign a confidence level to each prediction." (Predicting

In the case when the demands of this definition are not fulfilled the term "earthquake

We use such definition instead of another one when the term "prediction" means a deterministic prognosis as it was formulated in (Operational Earthquake Forecasting: State of Knowledge and Guidelines for Utilization. International Commission on Earthquake Forecasting for Civil Protection, http://www.protezionecivile.gov.it/cms/attach/ ex\_sum\_finale\_eng1.pdf). We suggest that deterministic prognosis is impossible now and hardly will be possible even in future, thus such use of the term "prediction" seems to be

**4.1 Case 1 - Diagnostics of a dangerous period before the 1994 Mw 8.3 Shikotan** 

1986, 1990) which provided a suitable procedure for prediction of large earthquakes.

functions used are following (Keilis-Borok & Kossobokov, 1986, 1990):

et al., 1980)) describes an increase in seismic activity;

B – describes the bursts of aftershocks.

within the circular areas with a fixed radius.

follows the concept of K. Mogi (1985).

L – describes deviation of N from the long-term trend value; Z – describes a linear concentration of earthquake sources;

The region under study in case 1 includes the Kuril Islands zone and the area to the east of Hokkaido Island. In 1992 we have prepared in computing form the earthquake M ≥ 4.0 catalog of the Kuril-Okhotsk region for the period 1962-1990. It was formed on the basis of yearly publications (The earthquakes in USSR…, 1964-1991). During the next years the catalog was updated by the Operative catalog data of the Sakhalin Branch of Geophysical Survey of the RAS. We used this regional catalog for testing the M8 algorithm (Keilis-Borok & Kossobokov,

The intermediate-term earthquake prediction technique, named M8 algorithm, is based on an assumption that a number of functions, defined for a particular earthquake sequence, become extremely large in values, within several months prior to a major shock. The

N – cumulative number of main shocks (aftershocks are excluded according to (Keilis-Borok

All functions, except the last one, were calculated twice: for a standard variant of small statistics (10 events or less per year) and for a standard variant of large statistics (20 events or more per year); where the numbers of events change by choice of threshold of magnitude taken into account. Two statistics are used for increasing robustness of results of prognosis. Values of these seven functions were used for adjusting the M8 algorithm, and then for diagnostics of Time of Increased Probability (TIP) for large earthquake (M ≥ 7.5) occurrence

Besides the method described above, we used a visualization technique to display spacetime distribution of seismicity to detect seismic gaps of the second kind. A gap of the second kind (seismic quiescence) refers here to a portion of a seismic area of low seismic activity with no observed earthquakes with М≥6.0 for a period of several years. This approach

aggregation in 0.1-bins is small, it influences the *b*-values weakly and uniformly; therefore it was not taken into account. The maximum likelihood method (3) gives a suitable *b*-value estimation for a number of events exceeding 50. Having this in mind the groups consisting of 50 subsequent events were used in *b*-value determination. The data points in Fig. 3 reflect the *b*-values obtained for such groups with step 25 events; thus the data points are independent of those next to the adjacent ones.

As it can be seen in Fig. 3, there is an evident tendency of decrease in *b*-values in the time vicinity of the generalized main shock; and this decrease increases strongly with approaching the moment of the main shock. In the foreshock sequence the noticeable decrease begins about one hundred days before the main shock. In the aftershock sequence the sharp increase in *b*-values takes place during the first several days after the main shock. A slow increase in *b*-values takes place in the following 100 days. It is necessary to notice that the *b*-values appear to be increased in comparison with the background value in the time interval 10-100 days after the main shock occurrence. These features agree with a tendency of lowermost *b*-values in the very beginning of the aftershock sequences and with an increase of *b*-value in the further evolution of the aftershock sequences (Rodkin, 2008; Smirnov & Ponomarev, 2004). The similar tendency was found in the examination of acoustic emission data (Smirnov & Ponomarev, 2004). New findings consist in the stronger decrease than it was found before and in rather symmetrical character of this decrease for fore- and aftershock sequences. Note that the amplitude of the *b*-value decrease appears to be proportional to the logarithm of time remaining from the moment of the main shock. Such type of behavior is typical of critical processes.

Note however that the well known and widely used below effect of "seismic quiescence" was not found in the generalized vicinity of strong earthquake. It can be connected with anisotropic character of this type of precursor anomaly in relation to a strong earthquake epicenter that is mentioned in (Zavyalov, 2006). In this case this effect can be eliminated by summarizing data from vicinities of a large number of differently oriented strong earthquakes.

#### **4. Experience in earthquake prediction at the Sakhalin Island and surrounding areas**

Region under study includes the Sakhalin Island and the Kuril Islands arc. In a few cases the area of the Japan Islands was also taken into account. This territory belongs to the transitive zone between the Pacific and the Eurasian continent and includes the active island arc characterized by one of the highest levels of seismicity on the Earth. Because of variability in quality of available catalogs the methodology of prognosis is more or less different in every particular case of strong earthquake prognosis, which is described below.

To avoid misunderstanding and controversial interpretations, we follow below the definition of the term "earthquake prediction," which was formulated by the Panel on Earthquake Prediction with the US National Academy of Sciences (Allen et al., 1976):

"An earthquake prediction must specify the expected magnitude range, the geographical area within which it will occur, and the time interval within which it will happen with sufficient precision so that the ultimate success or failure of the prediction can readily be judged. Only by careful recording and analysis of failures as well as successes can the eventual success of the total effort be evaluated and future directions charted. Moreover, scientists should also assign a confidence level to each prediction." (Predicting earthquakes …, 1976).

In the case when the demands of this definition are not fulfilled the term "earthquake forecasting" is used.

We use such definition instead of another one when the term "prediction" means a deterministic prognosis as it was formulated in (Operational Earthquake Forecasting: State of Knowledge and Guidelines for Utilization. International Commission on Earthquake Forecasting for Civil Protection, http://www.protezionecivile.gov.it/cms/attach/ ex\_sum\_finale\_eng1.pdf). We suggest that deterministic prognosis is impossible now and hardly will be possible even in future, thus such use of the term "prediction" seems to be inefficient.
