T. Parsons, S. Toda, R. S. Stein, A. Barka and J. H. Dieterich,
Heightened odds of large earthquakes near Istanbul:
An interaction-based probability calculation
Tom Parsons, Shinji Toda, Ross S. Stein,
Aykut Barka, James H. Dieterich
We calculate the probability of strong shaking in Istanbulan urban center of 10 million peoplefrom the description of earthquakes on the North Anatolian fault system in the Marmara Sea during the past 500 years, and test the resulting catalog against the frequency of damage in Istanbul during the preceding millennium. Departing from current practice, we include the time-dependent effect of stress transferred by the 1999 moment magnitude M=7.4 Izmit earthquake to faults nearer to Istanbul. We find a 62±15% probability (one standard deviation) of strong shaking during the next 30 years and 32±12% during the next decade.
The 17 August 1999 M=7.4 Izmit and 12 November 1999 M=7.1 Düzce earthquakes killed 18,000 people, destroyed 15,400 buildings, and caused $10-25 billion in damage. But the Izmit event is only the most recent in a largely westward progression of seven large earthquakes along the North Anatolian fault since 1939. Just northwest of the region strongly shaken in 1999 lies Istanbul, a rapidly growing city that has been heavily damaged by earthquakes twelve times during the past 15 centuries. Here we calculate the probability of future earthquake shaking in Istanbul using new concepts of earthquake interaction, in which the long-term renewal of stress on faults is perturbed by transfer of stress from nearby events.
Stress triggering has been invoked to explain the 60-year sequence of earthquakes rupturing toward Istanbul [Toksoz et al ref](1-3), in which all but one event promoted the next (4). Although an earthquake drops the average stress on the fault that slipped, it also changes the stress elsewhere. The seismicity rate has been observed to rise in regions of stress increase and fall where the off-fault stress decreases (5,6). The M=7.4 Izmit earthquake, as well as most background seismicity (7), occurred where the failure stress is calculated to have increased 1-2 bars (0.1-0.2 MPa) by M³6.5 earthquakes since 1939 (Fig. 1A) (8). The Izmit event, in turn, increased the stress beyond the east end of the rupture by 1-2 bars, where the M=7.2 Düzce earthquake struck, and by 0.5-5.0 bars beyond the west end of the 17 August rupture, where a cluster of aftershocks occurred (Fig 1B). The correspondence seen here between calculated stress changes and the occurrence of large and small earthquakes, also reported in (9), strengthens the rationale for incorporating stress transfer into a seismic hazard assessment.
A probabilistic hazard analysis is no better than the earthquake catalog on which it is based. Global observations support an earthquake renewal process in which the probability of a future event grows as the time from the previous event increases . To calculate such a renewal probability, ideally one wants an earthquake catalog containing several large events on each fault to deduce earthquake magnitudes, the mean inter-event time of similar events, and the elapsed time since the last shock on each fault . Although such catalogs are rarely, if ever, available, Ambraseys and Finkel compiled a wealth of earthquake damage descriptions for events since AD 1500 in the Marmara Sea region (12-15). We assigned modified Mercalli intensities (MMI) to 200 damage descriptions (available online), and used the method of Bakun and Wentworth (16) to infer M and epicentral location from MMI through an empirical attenuation relation (17). We calibrated the relation against Marmara Sea events that have both intensity and instrumental data (18). Uncertainties in earthquake location were explicitly calculated from MMI inconsistencies and inadequacies.
Our catalog thus consists of nine M³7 earthquakes in the Marmara Sea region since 1500. For the six events that occurred before instrumental recording began in 1900, we selected the minimum magnitude falling within the 95% confidence bounds at locations associated with faults of sufficient length (19) to generate the event (Fig. 2). We estimated rupture lengths and the mean slip from empirical relations on M for continental strike-slip faults (20). The locations and geometry of faults in the Marmara Sea are under debate; we follow (19), which is based on seismic reflection profiles (Fig 2), and find four faults capable of producing strong shaking in Istanbul: the Yalova, Izmit, Prince's Islands, and central Marmara. Our catalog suggests two earthquakes on the Izmit fault (1719,1999), yielding an inter-event time of ~280 yr, and three on the Yalova fault (1509, 1719, 1894), permitting an estimate of ~190 yr (21). We infer one earthquake (May 1766) on the Prince's Islands fault and one (1509) on the central Marmara fault (Fig. 2). For these, we gauge inter-event times by dividing the seismic slip estimated from the catalog by the GPS-derived slip rate (22,23), yielding a ~210 yr inter-event time for the Prince's Islands fault and ~540 yr for the central Marmara fault. Thus at least two of the four faults are likely late in their earthquake cycles.
One way to validate the catalog magnitudes, locations, and segment inter-event times is to compare the relative abundance of small to large shocks through the b-value; another is to see if the seismic strain release from the catalog is consistent with the measured strain accumulation from GPS. The frequency-magnitude relation for our catalog yields b=1.1 by maximum likelihood (24), close to the global average (25). Over a sufficiently long time period, the moment release by earthquakes must balance the moment accumulation by elastic strain if aseismic creep is negligible. We compared the seismic slip rate represented by the catalog (23.5±8 mm/yr) to the observed slip rate measured by GPS across the North Anatolian fault system in the Marmara region (22±3 mm/yr) (quoted uncertainties are one standard deviation here and elsewhere) (26) (Fig. 3). For b~1, most of the moment is conferred by the largest shocks, so the consistency between GPS and catalog strain means that the size and location of the three M~7.6 events, as well as the number of smaller earthquakes, are plausible.
Perhaps the strongest test of the 500-yr catalog can be made by calculating the combined Poisson, or time-independent, probability predicted from the inter-event times for the three faults we regard as capable of producing MMIVIII shaking in Istanbul. This is the probability averaged over several earthquake cycles on each fault, and yields 29±15% in 30 yr. This can be compared to the Poisson probability calculated directly from the longer record of MMIVIII shaking in Istanbul during the preceding ~1000 years (AD447-1508). The older record gives the long-term frequency of shaking used in a Poisson calculation without knowledge of the earthquake locations. At least 8 earthquakes (27) caused severe damage in Istanbul between AD 447 and 1508 (12-14), translating into a 20±10% 30-yr probability, roughly comparable to that derived from our catalog. Thus the fault inter-event times estimated from the 500-yr catalog are consistent with the independent record of shaking in Istanbul during the preceding millennium.
We combined earthquake renewal and stress transfer into the probability calculation on the basis that faults with increased stress will fail sooner than unperturbed faults. Because two out of the three faults within 50 km of Istanbul are interpreted to be late in their earthquake cycles, the renewal probability is higher than the Poisson probability. Additionally, the permanent probability gain caused by stress increase is amplified by a transient gain that decays with time. The transient gain is an effect of rate- and state-dependent friction (28-30), which describes behavior seen in laboratory experiments and in natural seismic phenomena, such as earthquake sequences, clustering, and the occurrence of aftershocks. We estimated the duration of the transient decay directly from the times between triggering and rupturing earthquakes on the North Anatolian fault (Fig. 4A). Because parameter assignments used in the calculation are approximate, we perform a Monte Carlo simulation to explore the uncertainties (31). The resulting probability functions (Fig. 4B) exhibit a gradual rise as the mean time since the last shock on each fault grows, and a sharp jump in August 1999 followed by a decay. We find a 62±15% probability of strong shaking (MMIVIII; equivalent to a peak ground acceleration of 0.34-0.65g (32)) in greater Istanbul over the next 30 yr (May 2000-2030), 50±13% over the next 22 yr, and 32±12% over the next 10 yr (Table 1). Inclusion of renewal doubles the time-averaged probability; interaction further increases the probability by a factor of 1.3.
The twelve earthquakes that damaged Istanbul during the past 1500 yr attest to a significant hazard, and form the basis for a 30-yr Poisson, or time-averaged, probability of 15-25%. Because the major faults near Istanbul are likely late in their earthquake cycles (with no major shocks since 1894), the renewal probability climbs to 49±15%. We calculate that stress changes altered the rate of seismicity after the 1999 Izmit earthquake, promoting the M=7.2 Düzce shock and the Yalova cluster. Because the 1999 Izmit shock is calculated to have similarly increased stress on faults beneath the Marmara Sea, the interaction-based probability we advocate climbs still higher, to 62±15%.
References and Notes
1. I. Ketin, Bull. Min. Res. Explor. Inst. Turkey 72, 1 (1969).
2. A. A. Barka, Bull. Seismol. Soc. Am. 86, 1238 (1996).
3. M. N. Toksoz, A. F. Shakal, and A. J. Michael, Pageoph 117, 1258 (1979).
4. R. S. Stein, A. A. Barka, J. H. Dieterich, Geophys. J. Int. 128, 594 (1997).
5. R. A. Harris, J. Geophys. Res. 103, 24,347 (1998).
6. R. S. Stein, Nature 402, 605 (1999).
7. A. Ito, et al., Precise Distribution of Aftershocks of the Izmit Earthquake of August 17, 1999, Turkey, Eos Trans. 80, F662 (1999).
8. We calculate Coulomb failure stress (, ) where D is the change in shear stress on the receiver fault, m is the coefficient of friction, D is the change in normal stress, and Bk is Skempton's coefficient. Stress values are found by elastic dislocation in a half space (33); viscoelastic effects are neglected. We used a slip model for the Izmit earthquake developed from InSAR (radar satellite interferometry) (34); slip models of other earthquakes are from (35) and (4). Here the assumed friction coefficient is 0.2, as has been found for strike-slip faults with large cumulative slip (36,37). A 100-bar deviatoric tectonic stress with compression oriented N55°W (38) is used, under which optimally oriented right-lateral faults strike E-W except along the rupture surface
9. Hubert-Ferrari et al., Nature 404, 269 (2000).
10. Y. Ogata, J. Geophys. Res. 104, 17995 (1999).
11. We use two probability distributions. The Brownian passage time function of Matthews (39) takes the form , where m is the average repeat time and a is the coefficient of variation. The lognormal distribution (40) is also used for time-dependent calculations. No catalog is adequate to estimate the coefficient of variation of the inter-event time, so we use a conservative value of 0.5 (10, 41).
12. N. N. Ambraseys, C. F. Finkel, Terra Nova 3, 527 (1991).
13. N. N. Ambraseys, C. F. Finkel, Terra Nova 2, 167 (1990).
14. N. N. Ambraseys, C. F. Finkel, The Seismicity of Turkey and Adjacent Areas: A historical review, 1500-1800 (Muhittin Salih EREN, Istanbul, 1995).
15. C. F. Finkel, N. N. Ambraseys, The Marmara Sea earthquake of 10 July 1894 and its effect on historic buildings, Anatolia Moderna Yeni Anadolu VII (Bibliothèque de l'Institut Français d'Etudes Anatoliennes-Georges Dumézil, Paris, 1996), vol. 43.
16. W. H. Bakun, C. M. Wentworth, Bull. Seismol. Soc. Amer. 87, 1502 (1997)
17. The relation Mi=(MMIi+3.29+0.0206di)/1.68, where di is distance in km between intensity (MMI) observation and epicenter, was developed from 30 California shocks with both intensity and instrumental observations (16). The RMS fit to this relation is calculated for trial locations on a 5x5 km-spaced grid. We excluded felt reports (MMI<IV) and MMI>VIII observations were saturated to VIII because criteria for higher intensities involve observations other than building damage, and because for poorly constructed masonry, damage may be total at MMI=VIII.
18. We calibrated with the 1912 Ms=7.4 Saros-Marmara (360 intensities; (42)), 1963 Ms=6.4 Yalova (11 intensities; (43)), and 1999 M=7.4 Izmit earthquakes (185 intensities). For the 1912 and 1999 events, we randomly selected 50 sets of 25 intensities (the mean number for the historical shocks) to calculate epicentral and magnitude errors. This yields intensity centers within ±50 km (at 95% confid.) of the instrumental epicenters, and gives the correct M or Ms within ±0.3 magnitude units. Site corrections were not made because we find no tendency for epicenters to be pulled toward sedimentary sites, and because improvement was only found in (16) when detailed site geology was available
19. J. R. Parke, et al., Terra Nova , in press.
20. D. L. Wells, K. J. Coppersmith, Bull. Seismol. Soc. Amer. 84, 974 (1994).
21. The most recent event for the Yalova segment is 1894.6; Izmit segment, 1999.7; Ganos fault, 1912.7; Prince's Islands fault, 1766.7; central Marmara fault, 1509.8.
22. Working Group Calif. Earthquake Probabilities, Seismic hazards in southern California: Probable earthquakes, 1994-2014, Bull. Seismol. Soc. Amer. 85, 379-439 (1995).
23. C. Straub, H.-G. Kahle, C. Schindler, J. Geophys. Res. 102, 27,587 (1997).
24. K. Aki, Bull. Earthquake Res. Ins. 43, 237 (1965).
25. S. G. Wesnousky, Bull. Seismol. Soc. Amer. 89, 1131 (1999).
26. Includes earthquakes in 1509, 1556, 1719, 1754, 1766, 1855, 1857, 1863, 1877, 1894, 1953, and 1964 from (12-14) and (35).
27. AD 447, 478, 542, 557, 740, 869, 989, 1323.
28. J. Dieterich, J. Geophys. Res. 99, 2601 (1994).
29. J. H. Dieterich, B. Kilgore, Proc. Nat. Acad. of Sci. USA 93, 3787 (1996).
30. S. Toda, R. S. Stein, P. A. Reasenberg, J. H. Dieterich, J. Geophys. Res. 103, 24,543 (1998).
31. The transient change in expected earthquake rate R(t) after a stress change, can be related to the probability of an earthquake of a given size over the time interval Dt through a nonstationary Poisson process as . Here r is the background seismicity rate, Dt is the Coulomb stress change, a is a state/rate constitutive parameter, s the total normal stress, t is time, and ta is the transient decay duration. The transient probability change is superimposed on the permanent change, which results from an advance or delay in the expected time until failure caused by the stress change. Integrating for is the expected number of earthquakes in the interval Dt, N(t) yields , where rp is the expected rate of earthquakes associated with the permanent probability change. This rate can be determined by again applying a stationary Poisson probability expression as , where Pc is a conditional probability, which can be calculated using any distribution. In addition to the inter-event time and elapsed time on each fault, this technique requires values for the stress change on each fault (we use the average calculated stress change resolved on each fault surface), the transient decay (shown in Fig. 4A from data in (4)), and a stressing rate on each fault derived from the fault geometry and the observed strain rate (0.1 bar/yr) (4). We performed 1000 Monte Carlo trials to establish error bounds (44). The four parameters for the Monte Carlo simulations are drawn at random from a normal distribution with a shape factor of 0.25 about each estimate, except for the inter-event time for which the shape factor is 0.5. Alternating Monte Carlo trials were run with a Brownian passage time and lognormal distribution.
32. D. J. Wald, V. Quitoriano, T. H. Heaton, H. Kanamori, Earthquake Spectra 15, 557 (1999).
33. Y. Okada, Bull. Seismol. Soc. Am. 82, 1018 (1992)
34. T. J. Wright, P. C. England, E. J. Fielding, M. Haynes, B. E. Parsons, Eos Trans. 80, F671 (1999).
35. Nalbant, S. S., A. Hubert, G. C. P. King, J. Geophys. Res, 103, 24469 (1998)
36. P. A. Reasenberg, R. W. Simpson, Science 255, 1687 (1992).
37. T. Parsons, R. S. Stein, R. W. Simpson, P. A. Reasenberg, J. Geophys. Res. 104, 20,183 (1999).
38. C. Gürbüz, et al., Tectonophysics in press (2000).
39. M. V. Matthews, J. Geophys. Res. in press.
40. Working Group Calif. Earthquake Probabilities, U.S. Geol. Surv. Circular 1053 (1990).
41. Working Group Calif. Earthquake Probabilities, U. S. Geol. Surv. Open File Rep. 99-517 (1999).
42. N. N. Ambraseys, C. F. Finkel, Annales Geophysicae 5B, 701 (1987).
43. N. N. Ambraseys, Earthqu. Engin. and Structural Dynamics 17, 1 (1988).
44. J. C. Savage, Bull. Seismol. Soc. Am. 81, 862 (1991).
45. The combined probability expression P=1-(1-Pa)(1-Pb)(1-Pc) for faults a-c assumes independent sources of hazard, since we cannot include future interactions and for all but the most recent earthquakes, we cannot include past interactions.
46. G. C. P. King, R. S. Stein, J. Lin, Bull. Seismol. Soc. Amer. 84, 935 (1994).
47. We thank N. Ambraseys, T. Wright, E. Fielding, A. Ito, J. Parke, and C. Finkel for sharing their insights and preliminary results with us, W. Bakun for his code and his review, and J. C. Savage, and W. Thatcher, C. Straub, and S. Kriesch for incisive reviews. Support from SwissRe is gratefully acknowledged.
Fig. 1. (A) Stress change caused by earthquakes since 1900. Shown are the maximum Coulomb stress changes between 0 and 20 km depth on optimally-oriented vertical strike-slip faults (8,46). Seismicity recorded since installation of IZINET (1993-July 1999 (7)) has uniform coverage over the region shown. Calculated stress increases are associated with heightened seismicity rates and with the future epicenter of the 17 August 1999 Izmit earthquake; sites of decreased stress exhibit low seismicity. Before the 1999 event, two studies (4, 35) identified this site as having increased stress and thus hazard. (B) Izmit aftershocks are associated with stress increases caused by the main rupture (first 12 days from IZINET (7)), such as the Yalova cluster southeast of "Y," and the occurrence of the 12 November Düzce earthquake. Faults: Y-Yalova, P-Prince's Islands, M-Marmara, I-Izmit.
Fig. 2. Large historical earthquakes since 1500. Intensities (dots) were assigned from damage descriptions compiled by (12-15). Red dashed contours give the moment magnitude M needed to satisfy the observations for a given location (16), because the farther the epicenter is from the observations, the larger the M required to satisfy them. The confidence on location is governed by the relative intensities; magnitude is a function of absolute intensities. We assigned earthquakes to faults by minimizing M within the 95% confidence region (17,18). Faults labeled in lower panels: I-Izmit, Y-Yalova, P-Prince's Islands, M-Marmara, G-Ganos, NAF-North Anatolia fault.
Fig. 3. Seismic slip from the 500-yr-long catalog of Fig. 2 is summed in four transects across the North Anatolia fault system in the Marmara Sea. All known or estimated M7 sources are included (26). The mean seismic strain release rate balances the strain accumulation rate observed from GPS geodesy (23). Whether earthquakes in parentheses extend to a given transect is uncertain. '1766a' is May; '1766b' is August.
Fig. 4. A. Transient response to stress transfer. The thirteen M³6.8 North Anatolian earthquakes for which the stress at the future epicenter was increased by ³0.5 bars are plotted as a function of time. The earthquake rate decays as t -1 in a manner identical to aftershocks, as predicted by (28-31). B. Calculated probability of a M³7 earthquake (equivalent to MMIVIII shaking in greater Istanbul) as a function of time. The probability on each of three faults is summed (45). The large but decaying probability increase is caused by the 17 August 1999 Izmit earthquake. Background tracks the probability from earthquake renewal; 'interaction' includes renewal and stress transfer. Light blue curve gives the probability had the Izmit earthquake not occurred.