Seismic quiescence before the M 7, 1988, Spitak earthquake, Armenia

Abstract

A detailed analysis of the 35yr of seismicity between 1962 and 1997 using a gridding technique shows that the M7, Spitak earthquake of 1988 December 7 was preceded by a quiescence anomaly that started at approximately 1984±0.5, and lasted about 5±0.5yr, up to the main shock. This quiescence anomaly had a radius of about 20±3km, estimated from circular areas with 75 per cent rate decrease, centred at the point of maximum significance of the anomaly. The quiescence was clearly present in the aftershock volume during the 5yr before the 1988 main shock, but its statistically strongest expression was located 30km NW of the epicentre. This anomaly fulfills the association rules between precursory quiescence anomalies and main shocks, even for a tight definition, and is therefore proposed as a case of precursory quiescence. The largest value of the standard deviate Z, found by random selection of samples by gridding, was Z=14 for a time window of T_w=3yr, using a sample size of N=300 events. This makes this anomaly the strongest observed so far, and it is the first documented in an environment of continental collision. There are no false alarms exceeding in significance the precursor. The Armenian earthquake catalogue used for this study had 4600 earthquakes with M≥M_min=2.2 in the area bounded by 39.5° to 42°N/42.5° to 47°E. From the point of view of homogeneous reporting this is the best catalogue we have analysed so far. The limits of the data used and the density of the grid are dictated by the data, and have no influence on the results. The choice of free parameters does not influence the results significantly within the following limits: 100≤N≤500, 2≤T_w≤7, 2.2≤M_min≤2.8.

1 Introduction

Precursory seismic quiescence has been reported by many authors (e.g. Mogi 1969; Wyss & Burford 1985, 1987; Matsu’ura 1986; Wyss & Habermann 1988a; Taylor et al. 1991; Ogata 1992; Wiemer & Wyss 1994; Wyss, Shimazaki & Urabe 1996; Yoshida, Ito & Hosono 1996). However, it is not clear yet how common this phenomenon is, what its characteristics are, and how it is best measured quantitatively. The purpose of this paper is to document in detail the case history of precursory quiescence before the Spitak M7 earthquake by the gridding method (Wiemer 1996) and to verify the results obtained by the SEISMOLAP technique (Zschau 1996a). Quantitatively mapped case histories of seismic quiescence are needed to determine whether the duration of quiescence is a function of the magnitude and of the tectonic environment, whether the size of the quiescence is a function of magnitude, whether the source volume itself or the volume surrounding it shows quiescence, and whether false alarms exist.

The quiescence hypothesis, as formulated by Wyss & Habermann (1988b), postulates that the quiet volume overlaps the main-shock source volume. This hypothesis assumes that some main shocks are preceded by seismic quiescence, which is a decrease of mean seismicity rate, as compared to the preceding, declustered background rate in the same crustal volume. The rate decrease, which must be judged significant by some clearly defined standard, takes place within part, or all, of the source volume of the subsequent main shock. It extends up to the time of the main shock, or may be separated from it by a period shorter than that of the quiescence itself. The rate decrease takes place in all magnitude bands. The duration of quiescence before large earthquakes is expected to be 4.5±3yr, but since it seems to depend strongly on tectonic environment Arabasz & Wyss 1996) and perhaps on the loading rate, it is important to measure it in the case of continental collision, since no such example of seismic quiescence has been quantitatively documented yet.

This hypothesis is based on the assumption that the seismicity rate in large and small seismogenic volumes is constant, unless a change in process takes place. Major earthquakes redistribute stress locally and thus cause aftershock sequences as well as local seismic quiescence on occasion. Creep episodes or slow earthquakes on deep extensions or other parts of faults also redistribute stress and hence cause local increases and decreases of seismicity rate. Thus, we must document how constant the seismicity rate is in the study area in order to measure quiescence. We do this by classifying all deviations from the average rate by their significance. If we find that the most significant, or one of the most significant, deviations from a constant rate occurs near a main shock, in time and space, we may propose that a connection exists.

The current version of the SEISMOLAP hypothesis of precursory quiescence postulates that the preparation process to main shocks influences the crust to distances of several hundred kilometres, causing a measurable seismicity rate decrease. Using laws of decreasing amplitude with distance and increasing duration with main-shock magnitude, the local seismicity regime may be modelled according to this hypothesis (Zschau 1996b), once the correct relationships have been found empirically. For Spitak, Zschau (1996a) concluded in a preliminary, simpler analysis that a one-year long quiescence in a volume of about 100km radius around the epicentre preceded this earthquake.

Arakelian & Martirossian (1996) found that the three main shocks (M≥6.7) in and near Armenia during the last two decades occurred in crustal zones of larger than average seismicity rate. For the highly active crustal volume north of the 1988 Spitak main shock, they demonstrated a reduction of seismic activity during a two-year window, which they positioned to last from 1986 to 1988.

The Spitak M7 earthquake of 1988 December 7, located at 40.92°N/44.22°E, was a reverse slip on a plane dipping to the NE (Arefiev 1991; Graiser 1991; Borcherdt 1993; Balassanian 1995). The aftershocks occurred in the top 12–15km. They define a rupture length of 40km and a width of 20km, approximately (Arefiev 1991; Borcherdt 1993). The maximum displacement along the 27km surface rupture was 2m reverse and 1.8m right lateral Rogozhin & Philip 1991; Balassanian 1995).

2 Data And Method

The earthquake catalogue for Armenia and the surrounding areas (Fig. 1) was produced by the National Survey for Seismic Protection of the Republic of Armenia by integrating the data of several regional catalogues. Epicentres are calculated to the nearest 0.1°, and depths are poorly known, other than that the earthquakes are in the crust.

Heterogeneous reporting as a function of time can generate false alarms and impede reliable measurement of natural rate changes (Habermann 1987, 1991; Wyss 1991; Zuniga & Wyss 1995). Therefore, we searched the catalogue for reported rate changes that may be artificial, based on an analysis by the algorithm genas (Habermann 1983), which searches for all significant rate changes in a range of magnitude classes cumulated from above as well as from below. Times at which more than 10 magnitude classes exhibit rate changes are selected for investigation by the magnitude signature technique. These changes are then interpreted as being probably due to artificial or natural causes, based on criteria explained elsewhere (Habermann 1991; Zuniga & Wyss 1995). This method was tested and found effective in cases where changes in network operations were well documented (Habermann 1986; Wyss 1991).

In the Armenian catalogue the reporting rate changes detected were far milder than in any of the numerous American catalogues we have investigated. This is interpreted as being due to the Soviet system of methodological and consistent observatory operations, compared to an American system of frequent innovation in hardware, software and analysis methods. A similar conclusion was reached by Eneva, Habermann & Hamburger (1995), who found the earthquake catalogue of another Soviet Republic to be of high homogeneity. We concluded that it is not necessary to make changes to the Armenian catalogue and that it can be used for seismicity rate analysis for M≥2.2. This is confirmed by the observation that the cumulative numbers do not suggest significant changes of reporting as a function of time for M≥2.2 (Fig. 2b).

The magnitude of completeness, M_c, as a function of time and space is of primary concern for defining mean magnitude (b-values). We mapped M_c in detail, using the computer program zmap (Wiemer, Zuniga & Allmann 1995). The algorithm we used interprets the point of greatest derivative of the non-cumulative frequency–magnitude distribution as M_c. Once we obtained maps of M_c, we checked the reliability of the results by examining the frequency–magnitude plot for numerous large and small subareas of each map. In the case of Armenia, the maps of M_c produced by the algorithm were invariably reliable. Readers who wish to check our work are welcome to use zmap and the Armenian data set.

We found that in the central area covered by the catalogue (39.7° to 41.8°N/42.7° to 46.7°E), the magnitude of completeness is M_c=2.2 for the period 1962–1997 (Fig. 2a). The magnitudes are calculated from the Russian system of energy classes. For this reason magnitudes fall preferentially into bins of M=2.8 and 2.2. The steps thus introduced in the frequency–magnitude relationship can be noticed (Fig. 2a). Nevertheless, the frequency–magnitude curve is well described by a straight line down to M=2.2. North, and especially south, of the boundaries indicated above, M_c=2.8. Therefore, we selected the central part of the catalogue (Fig. 1) with the aforementioned limits for analysis. This focused our work on the source volume of the M7 Spitak 1988 main shock and the volumes surrounding it.

In order to verify that our results do not depend on the choice of M_c, we performed the entire analysis in great detail for M_c=2.8 and M_c=2.6 as well. The results of this laborious work and the M_c maps are not shown in this paper because it would quadruple the number of figures.

The gridding technique we used to map quiescence has been explained elsewhere Wiemer & Wyss 1994), so we only summarize it briefly. We evaluate changes in the seismicity rate as a function of time at every node of a dense grid with 0.05° spacing. At each node, the data analysed consist of the nearest N events (here N=100, 200, 300, 400, 500 in separate analyses), spanning the time covered by the catalogue. We place a window (T_w=1, 2, ..., 8yr) into the time-series at each node and move it through time, stepping forward by one sampling interval (two months). For each window position, we calculate the standard deviate Z, generating the function LTA(t) (e.g. Wiemer & Wyss 1994), which measures the significance of the difference between the mean seismicity rate within the window, R₂, and the background rate, R₁, here defined as the mean rate outside the window, but in the same volume. The Z-value is defined as

where S and n are the variances and numbers of samples, within and outside the window. The resulting matrix of Z-values constitutes the basis of our quantitative analysis of quiescence. The Z-values of the LTA(t) functions are plotted at the beginning of the window for which they are evaluated, because we wish to define the beginning of a possible significant seismicity rate change.

The radii of the crustal volumes sampled are variable and inversely related to the local seismicity rate because we keep N constant to satisfy assumptions for statistical evaluations. In the labels of the Z-maps we give r_max, the maximum radius allowed. Nodes in low-seismicity areas, where r(N)>r_max, are not coloured on the maps because the information on seismicity rate is not local enough for analysis. r_max is found along parts of the edges of the coloured areas. Within the strong seismicity northwest of the Spitak epicentre (Fig. 1) the radii are at their lowest at about r=10km, typically.

The Z-maps we construct show the relative strength of rate changes at every node at the time of the window position selected. In movies we can view all possible window positions, and a summary of all high values of Z (above a selected alarm level) is given in alarm-cube representations of the results (Wiemer 1996; Wyss et al. 1996; Wyss, Console & Murru 1997). Thus, we can find all false alarms that may be present in a data set. To verify the reliability of a rate change mapped in this manner, we examine the cumulative numbers of earthquakes as a function of time in that volume.

For judging the significance of Z-values calculated for this data set, we generate synthetic Z-value matrices and the distribution of maximum Z-values, Z_max, that result from random selection of samples. The Armenian catalogue is the basis for this analysis. We extract from it, at random, N events and calculate the LTA function for a given window length (e.g. T_w=5yr). Repeating this process 1000 times results in a matrix approximately equal in size to the one obtained in the quiescence analysis. From this synthetic matrix we extract the Z_max value, and then we repeat this procedure 1000 times. The significance of a Z-value obtained in the analysis of the real data is judged in comparison to this synthetic Z_max distribution (Wiemer 1996; Wyss et al. 1996).

Some parameters of the analysis are free, i.e. they are chosen by the analyst and may be optimized to enhance the signal sought. Other parameters are determined based on factors unrelated to the analysis and cannot be changed by the analyst. Parameters that are set by circumstances not controlled by the analyst are the limits of the data set. The beginning of the modern earthquake catalogue in the study area is January 1962 (1962.0). During the first decade, the reporting was approximately homogeneous down to M=2.2, in the central area. Recent improvements of the resolution were not significant enough to warrant a later starting time for analysis with a lower magnitude level. The maximum number of usable earthquakes is obtained by starting in 1962.0 with a cut-off of M≥2.2. The end of the data set is 1997.03, the most recent data available. Thus, these limits of the data are not free parameters.

The parameters of node separation, bin length for individual rate estimates and the time step for the window advance are also parameters that are dictated by the nature of the data. Although they are selected by the analyst, and can be varied, they do not significantly influence the results. We select these parameters such that a nearly continuous coverage in space and time results, without excessive computing time. Thus, the node separation was selected as 0.05°, which is half the resolution of the epicentres. The bin length was chosen as two months, resulting in an average number of two events per bin in the analysis using N=300. It follows that the smallest step possible for moving the time window is two months. Thus, these parameters lead to coverage that can be described as nearly continuous; a denser coverage in space and time would not be warranted. The matrix that is generated to map the quiescence by Z-values contains approximately 8.4×10⁵ values. The spatial grid has about 4000 nodes, and there are 210 positions of T_w at each node.

Parameters that are varied to enhance the signal (quiescence) are the number of events, N, used in the selection at each node, and the window length, T_w, used for comparison of the mean rate in a window with the background mean rate. For results to be accepted as meaningful, we demand that they do not depend significantly on the choice of N and T_w, within a reasonable range. The duration of the quiescence is one of the parameters we must determine to define the anomaly. It is equal to the value of T_w that maximizes the significance of the anomaly. Likewise, the spatial extent of the anomaly must be measured, which can be done by finding the value of N that optimizes the significance of the quiescence. In this work we varied N between 100 and 500 events and T_w between 1 and 8yr. Small N and short T_w afford good resolution for defining quiescence, but the resulting statistical values are weak. For medium-sized N and T_w the statistical values are stronger, as long as the windows in time and space do not significantly exceed the extent of the quiescence. For very large N and T_w the statistical significance decreases again because the windows include data from areas with normal seismicity rates, diluting the quiescence signal. The table and figures show the results for medium-sized windows only.

The presence of clusters in the catalogue poses the following problems. If they are not removed (Reasenberg 1985), or accounted for by modelling (Ogata 1992) using some form of the Omori law (Omori 1900), they can generate erroneous alarms. These occur when episodes of normal seismicity rates follow cluster activity, and are thus mistakenly judged as quiescence. However, the modelling or declustering may also be imperfect, such that artefacts may be introduced into the data. In this paper we deal with this problem by (a) declustering the catalogue, using the Reasenberg (1985) algorithm with the constants proposed by Arabasz & Hill (1996), and (b) requiring that a quiescence anomaly, if found in the declustered data, must also be visible in the original catalogue including clusters. The total catalogue for the study area contains 5300 earthquakes with M≥2.2, including clusters. The declustered catalogue, the basis of our analysis, contains 4600 events with M≥2.2.

3 Analysis

The distribution of seismicity in the study area is diffuse (Fig. 1). The fault along which the 1988 Spitak rupture occurred does not make itself known by enhanced seismicity rates. The volume with about 20km radius, located adjacent to the Spitak rupture and northwest of it, displays an unusually strong seismicity rate, even in the declustered catalogue.

The Z-maps shown in Fig. 3 present different views of the data to demonstrate the differences and similarities of the anomalous seismicity pattern for different choices of N and T_w. Three maps are based on the declustered data, and one (Fig. 3d) on the original catalogue. The numbers of earthquakes selected around each node are N=400, 500 and 200 in Figs 3(a), (b) and (c), respectively. The time windows, in which the mean rate is compared to the mean background rate, are 5, 4 and 3yr respectively, and the starting times for the windows are 1983.9, 1984.9 and 1984.1 in the three frames. In all four maps, volumes north of the main shock are anomalously quiet during the selected time windows. Thus, these maps suggest that there may be a quiescence anomaly associated with the 1988 Spitak main shock. However, at first blush, the spatial correlation appears not to be as close as required by the hypothesis of Wyss & Habermann (1988b) and the temporal correlation has to be checked as well.

The cumulative curves (Fig. 4) of the number of earthquakes as a function of time show seismicity rate changes at specific nodes. These nodes were positioned randomly by the gridding process, but they were selected for presentation because they are located in the anomalous areas mapped in Fig. 3. Fig. 4 presents the quiescence anomaly in a number of different views to demonstrate that the quiescence is present, regardless of the choice of N or T_w. Declustered data for nodes using 500, 400, 300, 200 and 100 events are shown in Figs 4(a), (b), (c), (d) and (e). Cumulative curves for the aftershock volume (as defined by one month of aftershocks) for the clustered and the declustered data set are shown as Figs 4(f) and (g), respectively. In addition, a data set from a node of the clustered catalogue with 300 events is shown in Fig. 4(h). All of the cumulative seismicity curves show a clear anomaly before the 1988 Spitak main shock.

The statistical functions, LTA, displayed in Fig. 4 are the Z-values obtained by a comparison of the mean rate within a sliding time window and the long-term mean, defined by all seismicity outside of that window, but in the same volume. The maximum values of these functions are attained within the 6 to 2yr before the Spitak earthquake, and the quiescence lasts up to the main shock. The maximum Z-values of these LTA functions are large in comparison to those of other case histories (Wyss 1997), suggesting a high significance of this quiescence anomaly, which is correlated in time with the 1988 December Spitak M7 earthquake.

The fact that the aftershock volume also displays a quiescence anomaly in the cumulative number display (Figs 4f and g) is important because it shows that the main-shock source volume overlaps the anomalous volume. Although the aftershock volume produced no earthquakes during about five years before the Spitak main shock, it is not coloured by the darkest red in the Z-maps of Fig. 3 because the seismicity is low. A selection of N>100 earthquakes at a node located along the rupture reach from the source volume into parts of the crust that did not experience quiescence. Thus, the aftershock volume is coloured blue, green and orange in the maps (Figs 3c, d and b, respectively), although almost complete quiescence existed there.

After the 1988 Spitak earthquake, the declustered seismicity rate in most volumes returns approximately to the pre-quiescence level (Fig. 4). This observation strengthens the temporal correlation of the quiescence anomaly and the main shock. The temporal and spatial correlation of the quiescence (mapped in Fig. 3 and characterized in Fig. 4) with the Spitak main shock may not have any significance, if similar quiescence anomalies occur frequently at locations in space and time where no main shocks exist. Therefore, we search the matrix of Z-values (generated by the LTA functions at all nodes) for high Z-values, which may approach, or exceed, the Z-values recorded by the anomaly before the Spitak earthquake. The display of the high Z-values that could be false alarms is made in the alarm cubes (Wiemer 1996; Wyss et al. 1996, 1997). In these 3-D figures (Fig. 5) the bottoms are the spatial coordinates of the study area and the vertical axis is time. Alarms are defined as instances of Z-values larger than the selected alarm level at any node at any time. Fig. 5 shows that the results are stable, regardless of the choice of N (ranging from 200 to 500) and T_w (ranging from 3 to 5yr). In all alarm cubes for T_w≥3yr, there is only one group of anomalously high reports: before and near the Spitak main shock. For T_w<3yr, an alarm group, starting in 1994.7 and located near 41.4°N/43.8°E, shows Z-values approximately equal to those before the Spitak earthquake.

The relative strength of the Spitak quiescence anomaly compared to the strongest candidate for a false alarm (a quiescence equal to or greater than that of the Spitak alarm, but not followed by a main shock) can be gauged from Fig. 6, where the numbers of alarm groups are plotted as a function of the alarm level, Z_alarm. An alarm group is defined as a group of neighbouring nodes at which an alarm is declared at the same time. Fig. 6 shows the results using N=500 in the top row. N decreases in each row with N=200 at the bottom. In each column the window length is constant. The results using the largest, T_w=5, are on the left; those using T_w=3 are on the right.

The alarm group with the highest Z-value is always that associated with the Spitak earthquake, and it occurs well above the second strongest alarm. The Spitak quiescence anomaly is strongest compared to potential false alarms if N and T_w are large (e.g. Figs 6a and d). In these cases a large gap separates the Spitak anomaly’s Z-value from that of the alarm with the second largest Z-value, and few alarms reach Z-values above 5.5. If N and T_w are small, the separation between the largest and second largest Z-value of alarm groups becomes smaller, and many alarms reach Z-values above 5.5 (e.g. Figs 6i and l).

Quiescence episodes of the high significance discussed here do not occur frequently in the data set. The fraction of the total space–time volume occupied by alarms with Z>8 is typically below 1 per cent.

The statistical significance corresponding to a maximum Z-value is estimated by generating synthetic distributions of Z_max for every combination of N and T_w we used. The examples (Fig. 7) show that in all cases the value of the observed Z_max before the Spitak 1988 main shock was larger than any Z_max calculated by the random process. This indicates that the quiescence before the Spitak earthquake was significant above the 99 per cent confidence level.

4 Discussion And Conclusions

The Armenian earthquake catalogue is of unusually high quality from the point of view of homogeneous reporting in the area of Fig. 1. From 1962 to the present (35yr), there are no major changes of magnitude scale that we can detect. The reporting is not only homogeneous, but complete down to M=2.2 since 1962 (Fig. 2a). This means that the length of the catalogue available for defining background seismicity, and for searching for false alarms, is longer than other catalogues we have used, by a factor of two to three. As a result, the significance of the quiescence anomalies we find in Armenia (Z>10) is the highest we have ever measured (Wyss 1997). Thus, our conclusion that the 1988 M7 Spitak earthquake was preceded by precursory quiescence is uniquely firm. The hypocentre uncertainties in Armenia are larger than in most regional and local catalogues. However, this is of no consequence for the analysis of seismic quiescence, except that we are unable to map the quiescence as a function of depth in crustal cross-sections.

In the declustered catalogue, we find a single quiescence anomaly that far surpasses all other rate decreases in significance (Figs 4, 5 and 6). The centre of this anomaly is located about 30km north of the 1988 epicentre of the M7 Spitak main shock (Fig. 3). However, the anomaly also reaches into the 1988 main-shock source volume if viewed in large volumes, and it is clearly present in the aftershock volume (Fig. 4f). In the aftershock volume, the seismicity in the 20yr before the 1988 main shock was low: 57 earthquakes occurred. Two of these occurred a few months before the main shock and 54 occurred between 1962 and 1983.4. During the 5.3yr following 1983.4 only a single event occurred. Although this is a pattern of clear quiescence, it has low statistical significance, due to the low seismicity rate in this volume. Based on the 54 earthquakes that occurred before the quiescence, one would expect R₁=10 earthquakes during the 5.3yr of the quiescence, but only one occurred. Assuming that the earthquakes are produced by a Poissonian process, one would calculate the probability, p, of observing R₂=1 earthquakes as

which calculates as p<0.001. This evidence alone would not be sufficient to claim that a unique precursory anomaly existed before the Spitak main shock, because other, more significant anomalies may exist in the data set, but it is sufficient to show that the statistically highly significant anomaly documented to the north of the 1988 epicentre reaches into the main-shock source volume.

We defined the spatial extent of the quiescence anomaly by the area in which a 75 per cent rate decrease occurred. This criterion Wiemer & Wyss 1994) is arbitrary. In this paper, we define the area of 75 per cent rate decrease as follows: we increase the radius of a circle, centred at the coordinates of the largest Z-value, until the rate decrease, which is 100 per cent for small radii, is 75 per cent. If more than one node has the same Z_max, we select the one that leads to the larger anomaly volume. This method has the advantage that it is clearly defined and can be repeated. However, it has the disadvantage that it does not map the quiescence volume optimally because, by insisting on a circle as the anomaly shape, we will start to sample volumes in which no anomaly occurs in some directions, while in others the radius is not large enough yet to reach the edge of the anomaly. As a consequence, most circles defining the quiescence anomaly do not include the Spitak 1988 source volume, although it participated in the quiescence (Fig. 4f). Depending on the parameters chosen, the radius of the 75 per cent circular quiescence anomaly before the Spitak earthquake is 20±3km (Table 1). The number of events in the volume with 75 per cent rate decrease is not equal to N in the same row in Table 1; the latter is merely used for defining the centre of the quiescence with 75 per cent rate decrease.

We define the duration of the Spitak quiescence precursor, T_Q, as the time from the first report of the largest Z-value for the LTA(t ) function, t_Q, to the main shock. Onset time and duration vary somewhat as a function of the parameters chosen (Table 1). The average duration time is approximately 5±0.5yr and the quiescence starts at about 1984±0.5.

The quiescence anomaly before the Spitak earthquake is robust (Fig. 6 and Table 1). It shows up strongly in all the analyses we conducted, regardless of window length (1≤T_w≤8yr) and volume size selected around each node (100≤N≤600). In addition to the independence of the pattern on parameter choice, we find that it is present in the declustered as well as in the complete catalogue (Figs 3d, 4g and 4h). The maximum Z-value found by the random grid sampling ranges from 8.4 to 14.0, depending on the choice of N and T_w (Table 1). This makes this anomaly the strongest reported so far (Wyss 1997).

The minimum magnitude chosen for the analysis does not affect the conclusions reached for 2.2≤M≤2.8. If M=3.0 is used, the centre of the anomaly moves into the volume of the epicentre, but its statistical significance decreases because there are not many earthquakes above magnitude three and the quiescence can no longer be mapped. We conclude that the seismicity in Armenia is too low to define quiescence anomalies using large events (M≥3) only.

The location of the quiescence anomaly we mapped agrees with the location found by Zschau (1996a) and by Arakelian & Martirossian (1996). These authors report quiescence at some distance to the north of the Spitak 1988 epicentre. The fact that three independent methods find the quiescence further demonstrates the robustness of the anomaly. However, the relationship between the quiet volume and the main-shock source volume is complicated. As Fig. 5(f) shows, the source volume exhibited a nearly 100 per cent rate decrease during the 5yr before the main shock, yet the statistically strongest expression of the quiescence is located about 30km N of the epicentre (Figs 3a, b and d). When mapping the quiescence with small sampling radii (N≤200), one can see that the southern part of the seismically very active area (41°–41.4°N/43.7°–44.1°E in Fig. 1) does not turn quiet (Fig. 3c). Thus, it seems that constant seismicity production exists in a volume with about 5km radius, located between the quiet source volume to the south, and the statistically strongest quiescence to the north (coordinates of the samples in Fig. 4).

This feature and the strong expression of precursory quiescence in the hanging wall block of the Spitak reverse fault is predicted by the model of Kato, Ohtake & Hirasawa (1997). Based on a friction law, these authors modelled the approach to failure on a megathrust plane in a subduction environment. They found that precursory slip occurs on the plane and causes precursory seismic quiescence, especially in some parts of the shallow crust above the thrust plane. Their model describes accurately the seismicity pattern we found in the case of the Spitak reverse fault.

The duration of the Spitak quiescence and its uniqueness were documented in detail in this paper for the first time. The cumulative number curves (Fig. 4) clearly show that this quiescence lasted about 5yr. The alarm cubes (Fig. 5) and the plots of number of alarm groups as a function of alarm level (Fig. 6) show that the Spitak quiescence was uniquely strong, surpassing in significance all other alarms in windows larger than 2yr.

The Spitak quiescence anomaly is most strongly developed in data sets with N=300, 400 and 500 and T_w=3–5yr (Fig. 6 and Table 1). However, it is also unique (no other data sample produces a Z-value equal to or larger than that associated with the Spitak quiescence anomaly) for most other values of N and T_w, as summarized in Table 1. When the window length selected is short (T_w≤2yr for N=300 and 400) or the volume small (T_w≤3yr for N=200), a second anomaly, located adjacent to and northwest of the Spitak anomaly and lasting up to the present, rates approximately equal to the Spitak quiescence. It is generally true that shorter anomalies with smaller dimensions are more readily found in all data sets for which quiescence has thus far been mapped (Wyss et al. 1996, 1997).

We conclude that the 1988 M7 Spitak earthquake was preceded by a quiescence anomaly that started at approximately 1984±0.5yr, and lasted for about 5±0.5yr, up to the 1988 Spitak main shock. This quiescence anomaly had a radius of about 20±3km, as defined by the volume with 75 per cent rate decrease. Although it included the source volume of the 1988 rupture, its strongest expression occurred about 30km north of the source volume. This quiescence anomaly fulfils the association rules, even by the tight definition of Wyss & Habermann (1988b).

The extent of the anomaly is larger than the source volume of the main shock (Fig. 3 and Table 1). This is similar to the observations in Japan (Mogi 1969; Taylor et al. 1991; Ogata 1992; Katsumata & Kasahara 1996; Takanami et al. 1996; Wyss et al. 1996; Yoshida et al. 1996), where quiescences before M≥6.5 earthquakes extend to at least one source dimension beyond the aftershock volume, and in Utah Arabasz & Wyss 1996), where quiescences before M5 main shocks occupy volumes with about 20km radii. The large area of the Spitak quiescence anomaly also conforms to the SEISMOLAP hypothesis (Zschau 1996a). However, this observation is different from that before earthquakes along the San Andreas Fault system Wyss & Burford 1987; Wyss & Habermann 1988a; Wiemer & Wyss 1994) and in Hawaii (Wyss, Klein and Johnston 1981; Wyss 1986), where precursory quiet volumes have dimensions less than or equal to the main-shock rupture length.

The duration of the quiescence before the Spitak earthquake is similar to that before large earthquakes in California, the Aleutians (Kisslinger 1988; Ogata 1992; Kisslinger & Kindel 1994) and Japan, but it is also about the same as that before smaller main shocks in Utah Arabasz & Wyss 1996). If we assume that the duration depends on magnitude, the normal faulting regime of Utah with its slow loading rate seems to have exceptionally long precursor times. All other areas have approximately the same durations for precursory quiescence.

The detection of spatially disjointed quiescences is expected with the model of Kato et al. (1997), and it has been found in several tectonic environments. The first case where precursory quiet volumes separated by volumes with normal seismicity rate were found was that of the Kalapana M7.2, where a link between major asperities and quiescence was proposed (Wyss et al. 1981). In the case of the Landers M7.3 earthquake three separate precursory quiet volumes were found; two along the Landers rupture itself, and the third in the aftershock volume of the M6.5 Big Bear earthquake that followed the Landers main shock on the same day Wiemer & Wyss 1994). In the Armenia data set the second largest earthquake during the 35yr of recording was an M5.8 event that occurred at the northern edge of the strongest expression of the Spitak quiescence and half a year after the start of it (1984 September). It may be that this second largest earthquake in Armenia and the quiescence distribution were related to the same underlying tectonic process. Also, the 1988 Spitak earthquake occurred in a region where three major faults meet in a triple junction, and where diffuse seismic activity attests to the complexity of the tectonic setting. This is another reason why a complicated connection between the most quiet volume and the source volume may be expected.

To date we do not have enough evidence that would allow us to rule out strain softening (Stuart 1979) or dilatancy hardening (Scholz, Sykes & Aggarwal 1973) as the mechanism responsible for quiescence. Indeed, it may be that both mechanisms are operative in the Earth’s crust under different conditions. In the Kalapana and Landers cases the evidence favoured strain softening (Wyss et al. 1981) and dilatancy hardening Wiemer & Wyss 1994), respectively. We can only say that one should not be surprised that seismic quiescence is observed, because it is accepted that creep events on faults and pore-pressure changes occur. If these phenomena occur, seismic quiescence must follow them. Given the generally accepted model of strength heterogeneity of the crust (asperities intermixed with weak fault segments), one should also expect quiescence to occur in disjointed volumes. To design a firmly based model for seismic quiescence we need many more detailed case histories of quiescence, as well as much additional evidence from other geophysical parameters on the details of state and properties of the crust within which quiescence is observed and not observed. Also there is a hint of increased seismicity rate from about 1970 until the quiescence—a feature which seemed to exist in previously analysed cases as well. If this feature could be established at high significance levels, it might have implications concerning the mechanism of quiescence. Unfortunately it is not highly significant in any of the cases we have examined.

The fact that the spatial extent of the precursory quiescence anomaly at Spitak was larger than the aftershock volume and exhibited its statistically strongest expression about 30km north of the Spitak rupture shows that it will be difficult to issue accurate predictions based on quiescence anomalies that may be observed in real time. Up to now, we are aware of only three cases in which earthquake predictions were issued relatively successfully based on observations of quiescence. These were reported and discussed by Ohtake, Matumoto & Latham (1977, 1981), Wyss & Burford (1985, 1987), Kisslinger (1988), Wyss (1991) and Kisslinger & Kindel (1994). We believe that more experience in a number of tectonic environments is necessary to have a better idea of what to expect when seismic quiescence is observed in real time, and that observations of anomalies of parameters other than quiescence will be necessary before reliable predictions can be made.

Acknowledgements

This work was supported by the GeoForschungsZentrum and in part by the Wadati Foundation at the Geophysical Institute, University of Alaska, Fairbanks. We thank S. Wiemer for help with calculating the Z_max distributions and for critical comments on the manuscript. We also thank M. Garces, J. Lahr and J. Freymuller for comments on the manuscript.

References