Deterministic and Stochastic Models of the Seismic and Volcanic Events in the Santorini Volcano

In the second model, the observation that Thera is a simple Poissonian volcano, in conjunction with the Bayes theorem, leads to a probabilistic estimation of the date of the next eruption onset time. The third model is a particular case of the first one where the occurrence time of the next mantle shock(s) is determined from the observation that the process of occurrence of large mantle shocks in the South Aegean is Poissonian. All models imply that the next Thera eruptive cycle is not likely to start before the end of the present century.

INTRODUCTION

The Aegean and surrounding regions are very complicated from a seismotectonic and geophysical point of view (e.g. Papadopoulos et al. 1986), and the same might be said for some characteristics of the Cenozoic magmatism such as its petrochemistry and space-time evolution (Papadopoulos 1989a).

However, there are certain geophysical features with well-defined properties. One of the most important is the South Aegean calc-alkaline active volcanic arc. The caldera volcano of Thera (Santorini) is the most significant in this arc and one of the most significant in the Eurasian-Melanesian belt.

Thera has repeatedly erupted in post-Minoan times, right up to the present century. Studies concerning its possible activity in the future are of special interest in understanding how the volcano works and in mitigating the volcanic hazard. Only a very limited number of studies on the volcanic hazard and risk in this area have been published so far (Papadopoulos 1985; Fritzalas and Papadopoulos 1988).

As far as I know this is the first paper to be presented on forecasting the time of the next Thera eruptions. For this purpose three models have been tried. The first is a deterministic approach based on the suggestion that the onset time of eruption cycles depends upon the occurrence time of large subcrustal shocks. In the second model, which is purely stochastic, Thera is a simple Poissonian volcano. The third model is a particular version of the first model where the determination of the occurrence time of future subcrustal shocks is based solely on probabilistic terms.

THE DATA

All the known post-Minoan eruptions are described by Georgalas (1962; see also in Simkin et al. 1981, for the Volcanic Explosivity Index). There is evidence, however, that the reporting of volcanic eruptions is complete only since about the 16th century (Papadopoulos 1988). The dates of the eruptions which have occurred since that time, along with their characteristics, are listed in Table 1. This material has been used in the present analysis.

For the time interval 1900-1985 the earthquake catalogue of Comninakis and Papazachos (1986) is the data source for intermediate depth earthquakes analysed in the third model. This catalogue seems to be complete in events of surface-wave magnitude M_s ≥ 6.5. From the Bulletin of the Seismological Institute of the National Observatory of Athens it appears that intermediate depth events of such a magnitude were not recorded in the time interval 1986-1988.

A DETERMINISTIC MODEL

The close spatial coincidence of active volcanoes and inclined deep seismic zones was postulated many years ago (Gutenberg and Richter 1954). Subsequently, a wide variety of specific associations between deep seismic processes and volcanic activity in convergent plate margins has been proposed. One of the most interesting is the association of particular volcanic eruptions and particular intermediate depth or deep shocks preceding eruptions in several islands arcs including Thera in the Hellenic arc (Blot 1980). However, it remains questionable whether such an association is a certain feature of island arcs for reasons explained by Papadopoulos (1986, 1987a).

One of the problems of investigating spatio-temporal relationships between tectonic earthquakes and volcanic eruptions in a specified region is the sort of volcanic events which must be considered. For example, volcanic events such as steam or phreatic eruptions, smoke, and emission of 'more gas than usual' are inadequate for investigations of this type. Obviously it is desirable to consider only significant volcanic events which reflect directly the occurrence of magmatic processes in the magmatic chambers or even deeper in the roots of the volcanoes. For these purposes, magmatic eruptions which lasted at least 15 days and had a Volcanic Explosivity Index equal to or larger than 1, have been defined as 'significant volcanic eruptions' (Papadopoulos 1987b).

In a related paper, Papadopoulos (1987a) analysed data of large (M_s ≥ 7) subcrustal earthquakes (h ≥ 75 km), which occurred from 1900 to 1980, and found that in at least seven subduction zones of the circum-Pacific region periods of frequent volcanic activity with significant eruptions occurring in one to seven volcanoes of a given arc are systematically preceded by seismic energy peaks associated with usually no distant shocks of M_s ≥ 7.3 taking place directly beneath the volcanic belt of the same arc at the depth range of 80-190 km. The time lag between one energy peak and the onset of the following volcanic activity ranges between 1 and 18 years, but in most cases it does not exceed 12 years. In the same paper it is shown that on the basis of the x²-test in a p x q contingency table the correlation between the time series of large mantle shocks and subsequent eruption cycles is significant at the 99% level. The present author (Papadopoulos 1986) observed that in the South Aegean the occurrence of large intermediate depth shocks in the periods 1846-1863, 1887-1911 and 1926 (July)-1935 were correspondingly followed by the three eruption cycles which occurred on Santorini in 1866-1870, 1925-1926 (January) and 1939-1941, and include the minor event of 1950. The time lag is of the order of 15 years which is compatible with that found for the Pacific arcs. From these observations the suggestion that fracturing associated with large subcrustal shocks can cause magma formation or vice-versa is discussed in detail in the papers by the author cited above.

The observation that the significant volcanic activity in at least some island arcs may depend on the deep seismicity could contribute towards the investigation of an effective pattern for the long-term forecasting of volcanic eruptions. In the South Aegean it is noticeable that since February 3 1950, the day of the last Thera activity, neither volcanic activity nor intermediate depth shocks with M_s ≥ 6.8 have been reported. Assuming that the model mentioned earlier is correct, then, a new cycle of normal eruptions in Thera must be expected some years after the next large intermediate depth earthquakes beneath the South Aegean.

On this assumption the problem of forecasting the time of occurrence of the next eruption is equivalent to the problem of forecasting the occurrence time of the next large subcrustal shock(s). In a recent related paper the present author (Papadopoulos 1989b) utilized all the existing data and knowledge on this problem and, on the basis of seismicity observations and statistical considerations, concluded that such shocks are long overdue and that the next event(s) must be expected to occur in certain earthquake nests of the Benioff zone in the next few years.

In the next section a purely probabilistic approach to the problem is attempted. This approach, along with the deterministic model just discussed, generates a mixed model for volcanic eruption forecasting in Thera.

A STOCHASTIC MODEL

     Theory:     Earthquake occurrences can be modelled as a point process, that is a discrete stochastic process describing the positions of the events in time. The prototype of all point processes is the Poisson process which is defined by the three basic properties (e.g. Lomnitz 1974):

     (a) Independence    -    the number of events in any time interval is independent of the number of events in any other non-overlapping interval,

     (b) Orderliness    -   the probability of more than one event occurring in the same time interval is asymptotically negligible,

     (c) Stationarity    -   the probability of one event occurring in a short time interval dt is λdt, where is a constant in time.

The probability density function of x, the number of events per unit time, is the well-known Poisson distribution

f(n) = exp (-λ) λ ^x/χ!             (1)

Equivalently, the distribution of the time intervals T between successive events is the negative exponential distribution

f(T) = λ exp (-λT)                (2)

The probability of observing x events in time, t, is

P(x) = exp (-λt) (λt) ^x / x!    (3)

The rejection or adoption of the Poisson or other specific stochastic models depends on the selected time unit for which the numbers of events are counted. For example, Kárnik et al. (1983) examined the time distribution of the aftershocks of the M_s = 6.5, July 20 1978 event in the area of Thessaloniki, Northern Greece, and concluded that the negative binomial distribution is rejected for a one day time unit. On the other hand, it is not rejected for a time unit equal to only one hour.

To ensure that the selection of the time unit is as objective as possible, the present author has proposed (Papadopoulos 1989c) that the time unit must be taken as equal to the mean return period, T_m, of events having magnitude equal to or larger than the lower magnitude threshold, m, in the sample. The Gutenberg-Richter (G-R) magnitude-frequency relation provides an estimate of T_m. An independent estimate of T_m for large crustal earthquakes may be supplied from geotectonic methods, provided that proper data are available.

One of the most commonly used probabilistic approaches for the earthquake prediction problem is the Bayesian one. Although this is not a procedure on which specific predictions or special precautions can be based, it may still be quite useful for engineering purposes in the sense of long-term earthquake prediction. On the basis of Bayes's theorem in conjunction with the Poisson process model, Ferraes (1985) showed that the posterior Bayesian conditional probability, P(T_r/M), of the event T_r, given that an earthquake of magnitude M has occurred, can be expressed as:

P(T_r/M) = λT_r exp (-λT_r) [1 - exp (-λT_r)]                      

             ^{________________________________}                              (4)

               k

               Σ λ T_j exp (-λ T_j) [1 - exp (-λT_j)]

               j=1

This formula can be used to estimate the real time arrival dates for future earthquakes of a given magnitude range. Assuming that the inter-arrival time, T, has an exponential prior distribution and that we have n events each with inter-arrival time T₁, T₂,... , T_n, then, the parameter of the prior probabilities P(T_j) = 1 - exp (-λT_j), and conditional probabilities (likelihood function), P(M/T_j) = λT_j exp (-λT_j), is given by

         _{_}

λ = 1 T                          (5)

with

      _n

T = Σ T_i / n                    (6)

i = 1

The theory outlined above can also be applied to describe in probabilistic terms the expected times of the future volcanic activity in a specific volcano or volcanic belt under the condition that we take into account only the onset time of each eruption cycle regardless of its total duration. It is known that several 'simple Poissonian volcanoes' have been recognized in several regions of the world (see short review in Scandone 1983).

     The model:      To apply the previous theory to the Thera eruption cycles, two problems must be solved first:

1. How many eruption cycles took place in the time interval for which the data are thought to be complete, that is since the beginning of the 16th century?
2. Is Thera a simple Poissonian volcano?

The occurrence times of the known eight Thera eruptions since the 16th century AD are shown in Table 1. From their time separation it is evident that the events of 1570 (or 1573), 1650, 1707-1711 and 1866-1870 constitute independent eruption cycles. However, the frequent events of the time interval 1925-1950 indicate that probably more than one of the listed events constitute a unique eruption cycle. On the basis of certain volcanological features, Papadopoulos (1986) suggested that the 1928 and 1950 eruptions were the final events of the 1925-1926 and 1939-1941 activities, respectively. This analysis implies that six eruption cycles occurred in Thera in the examined time interval of 489 years. Nevertheless, in order to obtain results as reliable as possible two alternative procedures have been performed; one for six eruption cycles and another for eight cycles.

To investigate whether the volcano is Poissonian or not we have to examine whether the distribution of the eruption time onset follows the Poisson model. The determination of the probability density function, however, is difficult because of (a) the very small number of events in the sample, and (b) the lack of a reliable way to estimate the mean return period of eruption cycles, that is an objective time unit for which the number of events have to be counted.

To get over these difficulties I have examined how the distribution of the time intervals between successive eruption cycles fits the negative exponential distribution according to equation (2). The examination showed that the data fit the relation (2) for λ₁= 0.012823 events/yr determined by the least-squares method (Fig.1). The x²-test indicates that the Poisson model is adopted at the 0.01 level. The last question we have to answer is: what is the most appropriate value of λ to be taken into consideration? If we simply take λ = 6/489 we have λ₂ = 0.01227 events/yr. Alternatively, the mean of the time intervals, T, is given by (6). Here we have T = 73.83 (± 53.19) years. Hence, from (5) we get λ₃ = 0.013545 events/yr. Later on, in the Bayesian approach to the problem we have to use the value of Τ found above. For this reason λ₃ is also used in the present analysis. Moreover, λ₃ does not deviate practically either from λ₁ or from λ₂.

Now we can say that the volcano of Thera is Poissonian, at least in the last five centuries, and that the probability of observing x events in time, t, is

P(x) = exp (-0.013545t) (0.013545t)^x/ x!                  (7)

1-P(0), which is equal to 1-exp (-0.013545t), gives the probability of observing at least one event in time t. Table 2 shows values of P(0), 1-P(0), P(1) and P(2) for several t values. From the point of view of forecasting future events in Thera, the most important result in Table 2 is that the probability of observing at least one event in the next 50 years is about 50%, while the probability is being about 74% and 93% for the next 100 and 200 years, respectively.

The Bayesian approach to the problem reveals the real time arrival years of the next eruption cycles. Table 3 supplies information on these years as well as the corresponding Bayesian probabilities and time intervals. Each year is determined by adding each time interval to the year of onset of the last event, which is 1939. It is obvious that the highest Bayesian probability corresponds to the year 2019, while a past date has been 'predicted'.

Repetition of the same procedure for eight eruption cycles established that Thera is a Poissonian volcano (Fig. 2). The data fit equation (2) at the 0.01 level for λ₁ = 0.0183653 events/yr determined also in the least-squares sense. For reasons explained earlier the value λ₃ = 1 / Τ = 1 / 54.36 (±50.65) yrs = 0.018396 events/yr has been adopted instead of λ₁ and λ₂= 8/489 = 0.01636 events/yr.

After these modifications the probability of observing x events in time, t, is

P(x) =exp (-0.018396t) (0.018396t)^x/ x!                     (8)

Table 2 gives values of P(0), 1-P(0), P(1), and P(2) for several t values. Comparison with the corresponding values found in the case of six eruption cycles shows that in this case 1-P(0) is lower for the same time interval which means that the assumption of eight eruption cycles model increases the probability of observing at least one eruption in a given time interval in the future. This probability is about 60%, 84%, and 97% for the next 50, 100, and 200 years, respectively.

For the Bayesian approach, each of the seven time intervals is added to the year of onset of the last event, that is to 1950. Table 4 gives the real time arrival years of the next eruption cycle.

A MIXED MODEL

Here the first model, the deterministic one already discussed, is combined with a probabilistic approach to the expected seismic activity at intermediate depths beneath the South Aegean.

The approach is based on the theory outlined in the previous section. The statistical sample examined is composed by 15 intermediate depth earthquakes (h ≥ 70 km) of M_s = 6.6-8.0 which took place in the South Aegean area during the present century, that is in the time interval 1900-1988. The first and last of them occurred on August 11, 1903 and March 31, 1965, respectively.

The G-R relation shows that the mean return period of intermediate depth shocks of M_s ≥ 6.6 in the South Aegean is T_m = 8 years (Papadopoulos 1989b). This is the time unit selected to test whether the earthquake process is Poissonian. However, the number of events is too small to obtain a reliable probability density function of the number of events per unit time. To overcome the problem by extending the sample, a technique of largely overlapping time intervals, used by Papadopoulos and Voidomatis (1987) in similar problems, has been applied. Each time interval is equal to 8 years while the time shift is 1 year. The first 8-year interval is 1900-1907, the second is 1901-1908 and the last one is 1981-1988. A total of 82 intervals have been examined. From the total number of 15 events, 14 were counted 8 times and one event, that of 1903, was counted only 4 times. This means that 116 events have been counted and that each one of the real 15 events has been counted 7.73 times as an average. The mean rate of real events is λ₁ = λ₂/7.73 where λ₂ = 116/82 = 1.42 events/yr, that is λ₁ = 0.18 events/yr. Considering the mean number of events as the simple ratio of 15/89 we get λ₃ = 0.17 events/yr.

Observed and theoretical (expected) frequency distributions of the number, x, of the events analysed are shown in Fig. 3. According to the x²-test the distribution fits the exponential curve at the 95% confidence level which means that the process is being accepted as a Poissonian one:

P(x) = (1.42t)^x exp(-1.42t) / x!                             (9)

However, the probability of observing a real number of x events within time t is given by

P(x) = (0.18t)^x exp(-0.18t) / x!                              (10)

Table 5 gives the probabilities of observing 0, 1, 2, and at least one event in several time intervals. It is obvious that the probability of observing at least one event of M_s ≥ 6.6 in the next 25 years or so is about 99%. This result, in conjunction with the first model, leads to the conclusion that it is very probable to observe at least one eruption in Thera within the next 40 years approximately, provided that the first model is correct. This is consistent with the Bayesian elaboration of the volcanological data which implies that the years 2019 and 2030 have the highest Bayesian probabilities to signify the onset time of the next eruption depending on the number of past eruption cycles accepted.

In a Bayesian probabilistic prediction of intermediate shocks of M_s = 7.0-8.0 in the South Aegean, Papadopoulos (1987c) estimated as a rule past dates. This implies that if the first model is practically valid for earthquakes of M_s ≥ 7, then, the time interval extended up to next eruption must be narrowed down.

CONCLUSIONS AND DISCUSSION

The volcano of Thera was very active between 1925 and 1950. Since February 1950, however, it has been completely inactive leading thus to a lack of geophysical, geochemical and other observations useful in studying the evolution of the volcano and forecasting its future eruptions. For this reason, the problem has been approached by using deterministic and stochastic models based on the past history of the volcano.

The first model, the deterministic one, predicts that the onset time of the next eruptive cycle depends on the time of seismic energy peaks associated with large intermediate depth shocks. As shocks of this type are long overdue, we expect that such events may take place in the next few years. The time lag between seismic and volcanic activity is of the order of 15 years. The next eruptive cycle, therefore, must not be expected before the end of the present century.

According to the second model, which is purely stochastic, Thera is a simple Poissonian volcano, and the probability of observing at least one eruption increases with time. For the next 10, 25, and 50 years the corresponding probabilities are about 0.13, 0.29, and 0.49. A Bayesian approach indicates that no event is likely to occur before 1996. Due to incomplete data this analysis does not incorporate large eruptions which occurred before AD 1500. Geological records indicate that before the Minoan eruption another major volcanic event took place at 100,000 years BP (Heiken and McCoy 1984). The problem of incompleteness does not allow us to decide whether the process is Poissonian or time-dependent by taking into account the large magnitude events. Thus, we can say that the Poissonian model seems to be valid at least for eruptions of VEI = 2-3 (see Table 1).

From a data set of intermediate depth shocks of M_s ≥ 6.6 covering the post-1900 period it emerges that the earthquake process is Poissonian and that the probability of a new large shock in the next 25 years is about 99%. This means that a new eruptive cycle is very probable in about the next 15-40 years, provided that the first model is correct.

Previous estimations supply information on the expected future eruptive behaviour of Thera volcano. This picture is complemented by a qualitative volcanic risk determination which showed that the seismic component of the risk is the most important one not only because of the earthquake phenomena in the area but also because of the serious geotechnical problems in several residential zones mainly along the caldera edge (Fritzalas and Papadopoulos 1988). After these developments I feel that future research and administrative activities directed toward volcanic hazard mitigation must be mainly focused on (1) monitoring of the volcano, (2) elaboration of counter-measures, and (3) detailed study of the neotectonic and geotechnical conditions as potential factors of amplification of the seismic motion. The last is an urgent task because (a) the tourist activity in the island is very rapidly developing, and (b) the seismic component of the risk is associated not only with shocks or tremors of volcanic origin but also with tectonic earthquakes in the region of Cyclades, such as the large July 9, 1956 tectonic event of M_s = 7.5 in Amorgos which caused heavy damage, deaths and casualties on Thera.

Addendum

The problem of forecasting the time of occurrence of the next Thera eruption has been considered both in the author's paper, and also in the paper of Professor B.C. Papazachos. Comparison of the results must take into account that the two papers use slightly different data bases. Papazachos has taken the events of 1457, 1573, 1650, 1707, 1866 and 1925 as the largest ones in the post-AD 1500 interval. Thus, he found T = 93.8 years, sigma = 43.4 years and t = 64 years. However, the eruption in 1457 is questionable and, on the other hand, Papadopoulos additionally uses 1939-41 as a significant eruptive phase. These two changes imply that T = 73 years, sigma = 53 years and t = 50 years, which would change the probability of the occurrence of the next eruption in a given time interval determined in Papazachos' paper.

Another point relevant to both papers concerns the completeness of the data. Papazachos presented Fig. 1 which shows that the information available about volcanic eruptions in the Hellenic arc is complete after AD 1500. However, this diagram is just a repetition of that published by Papadopoulos (1988, Fig. 4), who was the first to supply evidence for the completeness of the data since the 16th century (see also in the section 'The data' above).

-------------------------------------------

----------------------------------------------