Examensarbete vid Institutionen för geovetenskaper
Degree Project at the Department of Earth Sciences
ISSN 1650-6553 Nr 393
Defining an Earthquake Intensity Based Method for a Rapid Earthquake Classification System
Definiera en intensitets-baserad metod för snabb klassificering av jordbävningar och förutsägelse av skador
Erik Bäckman
INSTITUTIONEN FÖR
Examensarbete vid Institutionen för geovetenskaper
Degree Project at the Department of Earth Sciences
ISSN 1650-6553 Nr 393
Defining an Earthquake Intensity Based Method for a Rapid Earthquake Classification System
Definiera en intensitets-baserad metod för snabb klassificering av jordbävningar och förutsägelse av skador
Erik Bäckman
ISSN 1650-6553
Copyright © Erik Bäckman
Published at Department of Earth Sciences, Uppsala University (www.geo.uu.se), Uppsala, 2017
Abstract
Defining an Earthquake Intensity Based Method for a Rapid Earthquake Classification System Erik B¨ackman
Ground motions caused by earthquakes may be strong enough to cause destruction of infrastructure and possibly casualties. If such past destructive earthquakes are analysed, the gained information could be used to develop earthquake warning systems that predicts and possibly reduce the damage potential of further earthquakes. The Swedish National Seismic Network (SNSN) runs an automated early warn- ing system that attempts to predict the damage of an earthquake that just got recorded, and forward the predictions to relevant government agencies. The predictions are based on, e.g. earthquake magnitude, source depth and an estimate of the size of affected human population. The purpose of this thesis is to introduce an additional parameter: earthquake intensity, which is a measure of the intensity with which the ground shakes. Based on this, a new earthquake hazard scheme, the Intensity Based Earthquake Classification (IBEC) scheme, is created. This scheme suggests alternate methods, relative to SNSN, of how earthquake classifications can be made. These methods will use an intensity database established by modelling scenario earthquakes in the open-source software ShakeMap by the U.S. Geological Survey.
The database consists of scenarios on the intervals: 4.0 ≤ Mw≤ 9.0 and 10 ≤ depth ≤ 150 kilometre, and covers the whole intensity scale, Modified Mercalli Intensity, 1.0 ≤ Imm≤ 10.0. The IBEC classi- fication scheme also enabled the creation of the ’Population-to-Area’ criterion. It improves prediction of earthquakes that struck isolated cities, located in e.g. valleys in large mountainous areas and deserts.
Even though such earthquakes are relatively uncommon, once they occur, they may cause great damage as many cities in such regions around the world often are less developed regarding resistance to ground motions.
Keywords: Ground motions, earthquake intensity, earthquake warning systems, damage prediction, classification
Degree Project E in Geophysics, 1GE031, 45 credits Supervisors: Bj¨orn Lund and Peter Schmidt
Department of Earth Sciences, Uppsala University, Villav¨agen 16, SE-752 36 Uppsala (www.geo.uu.se) ISSN 1650-6553, Examensarbete vid Institutionen f¨or geovetenskaper, No. 393, 2017
The whole document is available at www.diva-portal.org
Popul¨arvetenskaplig sammanfattning
Definiera en intensitets-baserad metod f¨or snabb klassificering av jordb¨avningar och f¨oruts¨agelse av skador
Erik B¨ackman
Markr¨orelser orsakade av jordb¨avningar kan va starka nog att skada v˚ar infrastruktur och orsaka d¨odsof- fer. Genom att analysera forna destruktiva jordb¨avningar och utveckla program som f¨ors¨oker att f¨oruts¨aga deras inverkan s˚a kan den potentiella skada minskas. Svenska Nationella Seismiska N¨atet (SNSN) driver ett automatiserat tidigt varningssystem som f¨ors¨oker f¨oruts¨aga skadorna som f¨oljer en jordb¨avning som precis spelats in, och vidarebefodra denna information till relevanta myndigheter. F¨oruts¨agelserna
¨ar baserade p˚a, t.ex. jordb¨avnings-magnitud och djup samt uppskattning av m¨ansklig population i det p˚averkade omr˚adet. Syftet med denna avhandlingen ¨ar att introducera ytterligare en parameter:
jordb¨avnings-intensitet, som ¨ar ett m˚att av intensiteten i markr¨orelserna. Baserat p˚a detta skapas ett jordb¨avnings-schema kallat Intensity Based Earthquake Classification (IBEC). Detta schema f¨oresl˚ar al- ternativa metoder, relativt SNSN, f¨or hur jordb¨avnings-klassificering kan g¨oras. Dessa metoder anv¨ander sig av en intensitets-databas etablerad genom modellering av jordb¨avning-scenarios i open source- programmet ShakeMap, skapat av U.S. Geological Survey. Databasen best˚ar av scenarior ¨over inter- vallen 4.0 ≤ Mw≤ 9.0 och 10 ≤ djup ≤ 150 kilometer, vilka t¨acker hela intensitetsskalan, Modified Mercalli Intensity, 1.0 ≤ Imm≤ 10.0. IBECs klassificeringsschema har ¨aven m¨ojliggjort skapandet av
”Population-mot-Area”-kriteriet. Detta f¨orb¨attrar f¨oruts¨agelsen av jordb¨avningar som tr¨affar isolerade st¨ader, placerade i t.ex. dalg˚angar i stora bergskjedjor och ¨oknar. ¨Aven om denna typ av jordb¨avningar ¨ar relativt ovanliga s˚a orsakar dom ofta enorm skada d˚a s˚adana h¨ar st¨ader ofta ¨ar mindre utvecklade r¨orande byggnaders motst˚and mot markr¨orelser.
Nyckelord: Markr¨orelser, jordb¨avnings-intensitet, jordb¨avnings-varningssystem, skadef¨oruts¨agelse, klassificering
Examensarbete E i geofysik, 1GE031, 45 hp Handledare: Bj¨orn Lund och Peter Schmidt
Institutionen f¨or geovetenskaper, Uppsala universitet, Villav¨agen 16, 752 36 Uppsala (www.geo.uu.se) ISSN 1650-6553, Examensarbete vid Institutionen för geovetenskaper, Nr 393, 2017
Hela publikationen finns tillg¨anglig p˚a www.diva-portal.org
Table of Contents
1 Introduction 1
1.1 Classification schemes . . . 2
1.2 Hypothesis . . . 4
2 Theory 5 2.1 Wave propagation, site amplification and source effects . . . 5
2.2 Ground Motion Prediction Equation . . . 7
2.3 Ground-Motion/Intensity Conversion Equation . . . 7
3 Data and software 10 3.1 Earthquake databases and data filtering . . . 10
3.2 Software and programming . . . 12
3.2.1 ShakeMap . . . 12
3.2.2 Generic Mapping Tools . . . 13
3.2.3 Python . . . 14
4 Method 15 4.1 Intensity simulations using ShakeMap . . . 15
4.1.1 Scripts for running ShakeMap . . . 15
4.1.2 GMPEs . . . 16
4.1.3 GMICEs . . . 22
4.1.4 Comparison to USGS . . . 22
4.2 The SNSN classification scheme . . . 25
4.3 The Intensity Based Earthquake Classification (IBEC) scheme . . . 25
4.3.1 Earthquake intensity . . . 27
4.3.2 Near-shore earthquakes . . . 30
4.3.3 Building standard . . . 31
4.3.4 Population . . . 32
4.3.5 Case studies . . . 32
5 Results 34 5.1 Evaluation of the SNSN scheme . . . 34
5.2 Evaluation of the IBEC scheme . . . 38
5.3 Comparison of the IBEC and SNSN schemes . . . 42
6 Discussion 45 6.1 Incorrectly classified earthquakes . . . 45
6.2 Differences in search radius . . . 46
6.3 Fault geometry and uncertainties in ShakeMap . . . 46
6.4 Future studies . . . 47
7 Conclusions 49
8 Acknowledgements 50
9 References 51
Appendices 52
A How to: ShakeMap 53
B Intensity function 56
C Distance to contour 57
D Haversine formula 59
1 Introduction
The strength with which the ground shakes while an earthquake occurs can be anything from not felt to violent enough to destroy buildings, and in the worst case cause casualties. Thus by investigating which factors that control the strength of the ground motions (also known as earthquake intensity) helps us to design programs that can quantify the earthquake intensities in real-time as data is being recorded by the seismographs. Further on, earthquakes may also be classified based on how much damage they are predicted to cause, based on the quantifications. The classification process is partly dependent on geological and geophysical properties but also society related parameters. Such parameters are ’how many people live in the vicinity of the earthquake epicenter?’ and ’what are the building standards in the area?’. Our programs tries to gather all information needed to classify the earthquakes and predict their shakings so that, in a next step, the society and concerned authorities are given as much and reliable information as possible.
This thesis focuses on hazardous earthquakes and how those can be classified based on their intens- ities. It also explores how, and if, we can predict the damage being caused by an earthquake that just happened, or for one that might happen in the future, using related parameters. The Swedish National Seismic Network (SNSN) records both local, usually low magnitude events, and tele-seismic events of moment magnitudes Mw≥ 4.0. As one of its tasks, the SNSN analyses all events Mw≥ 5.0 that occur globally and evaluates their risk to harm humans and infrastructure. The risk evaluation is presented to a few Swedish government agencies such as the Swedish Civil Contingencies Agency and the Swedish Ministry of Foreign Affairs. The analysis process is completely automated and SNSN is continuously running projects that are further developing the process to gain better precision and speed. This thesis is one such project and its main objectives are:
• Further develop the process of earthquake classification and damage prediction.
• Write a program that approximates earthquake intensities by tabulating outputs from the ShakeMap software.
• Evaluate the result of earthquake damage analysis done by the SNSN over the past years.
• Further develop the pre-existing functions for estimating population densities and classification of near-shore earthquakes.
These objectives will be summed up in the herein created Intensity Based Earthquake Classification (IBEC) classification scheme. The scheme is mainly dependent on earthquake intensity, from which it
also estimates size of affected surface area and population density, as being two important factors for a reliable earthquake classification.
Some words about the terminology. ’Event’ will be used for an earthquake that will, already have or might occur. ’Scenario’ is an imaginary event created by the user by choosing spatial coordinates and magnitude to model an arbitrary event occurring somewhere on the Earth. ’Intensity’ is sometimes used independently, without specifying what it is an intensity of, but it is always referring to the intensity by which the ground shakes during an event. The unit of intensity herein is based on the Modified Mercalli Intensity (MMI) scale which is explained in more detail in a separate section of this report.
1.1 Classification schemes
For earthquake classification and damage prediction the SNSN has developed a scheme that every re- corded tele-seismic earthquake goes through. The classification scheme requires four input parameters:
earthquake magnitude, hypocentral depth and epicentral latitude and longitude. Based on these paramet- ers the scheme will further determine an estimate of the affected population density, names and coordin- ates of nearby cities, whether the earthquake occurred inland or in an ocean, and closest plate boundary type that could have caused the earthquake. Thereafter the earthquake will enter the actual classification process and dependent on the values given to every parameter the earthquake will gain a class of either A, Bor C. The A class is the highest class in the sense that the damage following an A-event is assumed to be much larger than for a C-event. In fact the C class is given to earthquakes with low damage potential.
Figure 1 shows the magnitude-depth relationship for approximately 20 000 events recorded by the SNSN and the colouring represents each event’s class given by the SNSN classification scheme. Note how class Aand B extend to deeper depths as magnitude increases and the step-like behaviour seen, especially in the B class, is due to depth dependent classification boundaries in the SNSN scheme. The plot order is C→ B → A so there is an overlap with the highest class always on top.
The herein created IBEC classification scheme is a modification of the SNSN scheme and is using a foundation provided by the latter. It is a modification in the sense that the IBEC scheme does its classifications based on earthquake intensities rather than earthquake magnitude and depth, as intensity is, by the author, thought to be a better estimator of earthquake damage. The IBEC scheme also estimates affected population densities and distinguishes between inland and near-shore earthquakes in slightly different ways. IBEC also brings in and uses each country’s Gross Domestic Product (GDP) per capita (year 2015) as an estimator of building standard where a building with a high standard is more resistant to ground motions reducing the risk of a collapse. The global distribution of GDP per capita is seen in Figure 2 and is divided on countries below and above a GDP per capita of 20 000 International dollar (Int $).
Figure 1. Magnitude-depth distribution for approximately 20 000 earthquakes recorded by the SNSN. The colors represent the classes given by the SNSN classification scheme: A = high damage, B = moderate damage and C = low potential damage.
Figure 2. Global distribution of GDP per capita. Each country is coloured either red or blue depending the country’s current (year 2015) GDP per capita: red ≥ 20 000 Int$ and blue < 20 000 Int$.
In the result section a comparison of the performances by the two schemes is made. The SNSN scheme defines their classes to be A = high damage, B = medium damage and C = low potential damage, and unlike IBEC, does not relate its classes to any predicted death tolls. Thereby will the comparison be according to the definitions set by the IBEC scheme: A ≥ 100, 1 ≤ B < 100, and C = 0 death tolls.
To temporarily convert the SNSN classes into death toll numbers might be slightly uncertain, but for simplicity the comparison is preferably kept numerical.
1.2 Hypothesis
The earthquake data that are processed and analysed within this project are contained in two sets of real data. See Section 3.1 for a detailed review. The first data set will be used for evaluating the result of the past years earthquake analysis done by the SNSN. The second set contains data about the number of casualties caused by earthquakes over the past twelve years. This data set will be used when further developing the earthquake classification and damage prediction process as it gives information about the true aftermath of a particular earthquake. This casualty-data set, lets call it ’death toll data’, will be subdivided into three groups: A, B and C. Which group an earthquake belongs to depends on how many casualties that were caused by the earthquake and the distribution, according to IBEC, is A ≥ 100, 1 ≤ B < 100 and C = 0 casualties. These groups will be plotted versus different parameters in various ways to detect possible relationships. Then, an attempt of reproducing these subdivided groups and plots will be made by using, for instance, the intensity and estimated population and without using the actual death toll data. If a reproduction is possible, using any parameter other than death toll, then what is left is a model that can approximately predict the damage following an earthquake.
The methodology used for solving this problem is influenced by inversion theory. The death toll data in Section 3.1 is our observed data, ¯dobsand we seek a model, ¯m, which minimizes ¯dobs− ¯dpre= ¯dobs−G ¯m where ¯dpreis the predicted data and G is the data kernel containing all parameters that we can apply and think reflects what is seen in ¯dobs. What the final ¯mactually will consist of is a set of if-statements in Python that every incoming event will go through and depending on the paths taken by the earthquake it will be classified as an A, B or C earthquake where each class refers to the casualty intervals given earlier.
The classification scheme requires following parameters: earthquake magnitude, hypocentral depth and epicenter latitude and longitude.
2 Theory
2.1 Wave propagation, site amplification and source effects
A rupture along a fault plane will cause radiation of seismic energy that travel outwards from the fault as pressure and shear waves (P- and S-waves, respectively) and their joint name is body waves. Portions of body wave energy will convert into surface wave energy within the Earth’s crust and travel as Rayleigh- and Love-waves (R- and L-waves, respectively) with a propagation direction parallel to the surface.
Despite the fact that all four wave types cause ground motions, surface waves have the largest motions if it is measured and compared at the surface. Ground motions in general will decrease in amplitude, as source distance increases, through anelastic attenuation and geometrical spreading. Anelastic attenuation absorbs energy due to inhomogeneities within a material while geometrical spreading causes a decrease in energy density as the wave energy is gradually spread over larger volumes as the wave propagates.
One purpose of having seismographs recording wave amplitudes is that such amplitude data can be used to form Ground Motion Prediction Equations (GMPEs) which are used to predict the level of ground motions on the surface anywhere around a source. Thus, it is possible to predict the earthquake intensity for an event that just happened or for an event that may happen in the future, and GMPEs are therefore widely used in earthquake warning systems. GMPEs are introduced in more detail in Section 2.2 and how they are modelled using the ShakeMap software in Section 3.2.1. Working with GMPEs require knowledge, or at least an assumption, of wave velocities for the upper parts of the crust as ground shaking is amplified when waves travel through a low velocity material. Data bases of surface and near- surface geology where also S-wave velocities are known from e.g. bore hole measurements are generally hard to find. On the contrary, topographic data is globally well defined and relatively easily accessed.
Topographic variations could give an intuitive feeling about geomorphology where steep mountains can be considered as hard rocks whereas valleys would be filled with sediments. Wald and Allen (2007) found a mathematical relationship between the topographic gradient and S-wave velocity of the upper 30 meters of the crust, Vs30, by simply plotting observed Vs30data versus topographic gradient data. They had a significant scatter in their result but through application of statistical methods and referring to other studies concluding that S-wave velocity increases as grain size and void spacing increases, Wald and Allen (2007) claims that the scatter is within the acceptance margins. The relationship says that steeper topographic gradients give higher S-wave velocities and vice versa, which agrees with data from e.g. Mavko (2015) stating that a typical S-wave velocity for granites is 2500-3300 m/s and for dry sands is 100-500 m/s. Also, an increase in lithification of the sand result in a higher velocity.
In the real world, earthquake sources are plane sources with asperities rather than the common as- sumption of point sources. A drawback of assuming point sources is that, in the intensity modelling process, the greatest particle movements seems to appear in a point, at the epicenter. While the true fact is that the greatest movements occurs all along the whole rupturing fault plane and this model error increases with larger faults. Thereby knowledge about the dimensions and orientation of the rupturing plane increases the validity of the models. An example of when source dimensions are known has been modelled by the U.S. Geological Survey (USGS) and is seen in Figure 3. Intensity is highest in the vicinity around the surface trace of the fault and decreases with increased distance to the trace. Known dimensions and fault orientation also helps understanding the source effects, like directivity. Directivity occurs in the direction in which the rupture propagates, i.e. the source is moving along a line with time.
Imagine a source that is propagating towards a person. The energy radiating from the far-end may con- structively interfere with the energy radiating from the near-end and the person will experience a greater amplitude than if the same source was propagating away from the person.
Figure 3. Earthquake intensity modelling when source dimensions are known. The intensity is greatest along slightly curved line which is the surface trace of the fault. The star denotes epicenter. This ShakeMap is made by the USGS.
2.2 Ground Motion Prediction Equation
Ground Motion Prediction Equations (GMPEs) are empirical relationships used to explain the attenu- ation of earthquake ground motions as distance to epicenter increases. Thereby it is possible to predict how strong ground shaking one can expect for a given magnitude and depth at different locations around the epicenter. This is especially useful for areas with sparse station network coverage or if a poor na- tional economical situation prevents installation of any station. GMPEs are generally found through the application of regression analysis methods, i.e. finding a function that best fit the observed data. Such data are peak ground acceleration (PGA), peak ground velocity (PGV) and peak ground displacement (PGD) data which are measured at the seismographs, or rather accelographs as regular seismographs tend to saturate if they are located too close to an epicenter. PGA, PGV and PGD are sometimes lumped together in the term peak ground motions (PGM).
The attenuation relationships vary depending on for example source type (compressional, extensional or transform sources), earthquake magnitude, source depth and geology. Youngs et al. (1997) found that, in a compressional region, events that occur at the compressional interface where the subducting and overriding plate meet is different from events that occur within the subducting plate, so called intraslab events. Intraslab events are typically high-angle normal-faulting events. They also found that ground motions related to subduction zone earthquakes are stronger than motions related to shallow crustal earthquakes, and this pattern increases with magnitude. Their attenuation relationship is a function of earthquake magnitude, distance to rupturing plane, hypocentral depth and source type, where they distinguish between interface and intraslab events, and regression coefficients.
GMPEs have been established by many different researchers and groups. Each GMPE has its own characteristics and is most suitable for one kind of tectonic region. Table 1 lists the GMPEs that are considered and modelled in this project (see following sections) and all of those are already implemented into the ShakeMap software. There is also a possibility to implement new GMPEs into the software for future use.
2.3 Ground-Motion/Intensity Conversion Equation
Ground-Motion/Intensity Conversion Equations (GMICEs) are used to convert peak ground motion data to a value on the Modified Mercalli Intensity (MMI) scale, 1 ≤ Imm≤ 10 (see Table 2). Peak ground motion can be either acceleration, velocity or displacement while intensity is a single value representing all three PGMs together. The MMI scale greatly simplifies the representation of the ground motions that follow an earthquake to the society. As MMI is an intensity scale, unlike regular magnitude scales, it describes the observable effects of an earthquake such as brittle deformation of concrete walls. These
Table 1. The Ground Motion Prediction Equations (GMPEs), that are used within this thesis, for modelling earthquake intensity. Magnitude and distance range are the ranges of real data used when the reference established the equation. Region is where the real data came from.
GMPE name Reference Magnitude range Distance range (km) Region
AkkarBoomer07 Akkar and Bommer
(2007) 5.0 ≤ M ≤ 7.6 5 - 100 Europe
ASB13 Akkar et al. (2013) 4.0 ≤ M ≤ 8.0 0 - 200+ Southern Europe and Middle East
MA2005 Motazedian and
Atkinson (2005) 3.0 ≤ M ≤ 8.0 2 - 500 Puerto Rico
Youngs97 Youngs et al. (1997) 5.2 ≤ M ≤ 8.0 0 - 300 Subduction
effects, or damages, are reported by the society as ’Did You Feel It?’ (DYFI) data and for damage prediction purposes such data may be hard to calculate with just mathematics. Instead by correlating observed and reported DYFI data with PGM data measured by the accelographs one may achieve an empirical relationship between MMI and PGM. Like GMPEs there are different conversion equations established in different regions and the relationship found by Wald et al. (1999) is used through out this project. They found that instrumental recordings of peak ground acceleration tend to saturate at high intensities and therefore use both peak acceleration and velocity for the regression analysis. Their final relationships show that intensity depends logarithmically on peak ground motions. Dewey et al.
(2002) found that DYFI data were consistent with PGM data being converted to intensity over the whole MMI scale apart from some scatter in the lower intensities. Section 4.1.3 will discuss the ShakeMap implemented GMICEs in more detail.
Table 2. The Modified Mercalli Intensity (MMI) scale which relates peak ground motions with observable effects using numerals. Descriptions from e.g. USGS (2015b).
Numeral Description Numeral Description
I Not felt. VI
Felt by all. Many frightened. Plaster
might break.
II Maybe felt in upper
floor buildings. VII
Damaged negligible in buildings with good design. Slight damage
in well-built buildings.
Considerable damage in poorly built buildings. Some chimneys broken.
III Felt outdoors as if a
truck passes by. VIII
Considerable damage in well-built buildings. Great damage to poorly built
buildings, partly collapse.
IV Walls make cracking
sound. IX
Great damage in well-built buildings.
Disaster in poorly built buildings.
V
Many people awakened. Some windows broken.
X Total disaster.
3 Data and software
3.1 Earthquake databases and data filtering
The earthquake data presented herein come from two different databases. One set consists of approxim- ately 20 000 tele-seismic events occurring between year 2002 to 2015 gathered by the SNSN, whereas the other set of approximately 1000 events span magnitudes 1.0 ≤ Mw≤ 9.0 occurring between year 2000 to 2012, and comes from the ’Significant Earthquakes of the World’ database at the USGS (USGS (2015c)). The bigger set is partly used to evaluate the past years analysis done by the SNSN and to find appropriate boundary values used within the different functions that are written. The smaller data set covers earthquakes that have caused ground motions strong enough so that they have been felt and in most cases also caused casualties. If there have been casualties the amount is given and thus will be treated as observed data (see Section 1.2) when trying to further develop the earthquake classification process.
Events with known number of casualties will be divided into two groups. One group of 107 events where death toll ≥ 3 and one group of 318 events where death toll = 0. 156 events, where 1 ≤ death toll < 3, are excluded from the analysis because the casualty reason for these earthquakes is mostly heart attacks and those are not necessarily related to strong ground motions. Figure 4 show the geographical distribution of events, with death toll ≥ 3, gathered from the ’Significant Earthquakes of the World’
database. Common reasons for the casualties are rock and landslides initiated by ground motions, col- lapsing houses, heart attacks and other less common reasons. Figure 5 show the distribution of excluded earthquakes with death tolls of 1 or 2. Figure 6 show data from the same database but for events where death toll = 0.
Building codes specify the building standards required to withstand ground motions to a certain de- gree. A poor economical situation counters the possibility to apply building codes and data of GDP per capita will be used to approximate each country’s building standard. The data used comes from CIA (2015) and the method applied in IBEC is explained in Section 4.3.3.
Figure 4. Geographical distribution of 263 earthquakes with a death toll ≥ 3. Each marker is an earthquake of magnitude 3.8 ≤ Mw≤ 9.0 that occurred between the years 2000 and 2012. Marker size increases with death toll.
Figure 5. Geographical distribution of 263 earthquakes with a death toll between 1 and 2. Each marker is an earthquake of magnitude 2.2 ≤ Mw≤ 7.7 that occurred between the years 2000 and 2012. Marker size increases with death toll.
Figure 6. Geographical distribution of 763 earthquakes with a death toll = 0. Each marker is an earthquake of magnitude 1.0 ≤ Mw≤ 8.3 that occurred between the years 2000 and 2012.
3.2 Software and programming
3.2.1 ShakeMap
ShakeMap is an open-source program developed by Wald et al. (2005) and maintained by the USGS.
It runs together with several other programs and modules such as the Generic Mapping Tools (GMT).
For detailed instructions on how to install and configure the program see the ShakeMap manual. In Appendix A a short ’How To’ is presented of how scenarios are generated in ShakeMap.
ShakeMap is used to model various parameters such as earthquake intensity, peak ground motions and spectral acceleration. Spectral acceleration is used by earthquake engineers to study the oscillating behaviour of structures affected by ground motions. Within this project, ShakeMap is used for modelling the variability of earthquake intensity around an epicenter. Events studied with ShakeMap are referred to as scenarios, i.e. the user creates an arbitrary event by choosing spatial coordinates and earthquake magnitude. These scenarios are modelled using the most general settings possible to avoid tuning the software into any seismic locality as SNSN and IBEC analyses incoming global events. Such general settings assume a point source and a constant S-wave velocity of 686 m/s in the upper 30 m of the crust and no local site amplifications or source effects, regardless of earthquake magnitude and source depth.
If one wants a more detailed model of earthquake ground motions for a particular location, this can be achieved by adding a local site amplification grid file, the Vs30relationship by Wald and Allen (2007) discussed in Section 2.1. This relationship allows the user to estimate local near-surface geology and S- wave velocities anywhere using topographic data. The grid file is implemented into ShakeMap to more
accurately predict earthquake ground motions for the area of interest. Note that adding such grid files make every model locally unique. The well-known Tohoku-Oki earthquake that occurred on the 11th of March 2011 is modelled in Figure 7. These figures show the difference in resulting intensities between including a site amplification file and excluding the same file. Figure 7a models the event using a Vs30 layer and there are indications of intensity amplifications of as much as about 1.2, or an increase from I = 7.5 to I = 9.0, in the bay on the Japanese east coast. This increase indicates the potential damage to be heavy instead of moderate. It is seen that intensity amplifications exist where the topography is lower compared to its adjacent areas, such as at the coasts and in valleys. This pattern is to expect as the derivation of Vs30is based on topography variations (Wald and Allen (2007)).
(a) Site amplification grid included. (b) Site amplification grid excluded.
Figure 7. Comparing the resulting ground motions intensities between including and excluding a site amplification, Vs30, grid. Its seen that ground motions gets amplified in depressions filled with loose sediments such as valleys and bays. The event plotted is a re-creation of the Tohoku-Oki earthquake of magnitude 9.0 and a hypocentral depth of 29 km that struck Japan the 11th of March 2011.
3.2.2 Generic Mapping Tools
The Generic Mapping Tools (GMT) is an open-source software by Wessel and Smith (1995) used for manipulating and plotting Cartesian and geographical data. It is used globally and within for example earth science, physics, mathematics and biology. It is a program without user interface but it comes with a well-written program manual and online help. For the purpose of this project, GMT is used for plotting the distribution of different data on geographical maps. GMT is also used within ShakeMap as it makes contour maps from what is calculated by different equations, e.g. GMPEs and GMICEs.
3.2.3 Python
The programming language Python has been used throughout the whole project. It is a freeware and requires only a regular text editor and a command prompt to execute scripts. Python has many in-built functions and relies on row indentation rather than various characters, such as {, } and ;, to define blocks and row ends which makes Python a clean and easily read language. Also, many of the already existing scripts at SNSN are written in Python.
4 Method
4.1 Intensity simulations using ShakeMap
Referring to Section 3.2.1, ShakeMap have been used to create earthquake intensity matrices by tabulat- ing the results from modelling scenarios of different magnitudes and source depths. The following sub- sections deal with how the desired intensity data have been extracted from ShakeMap and also discusses the available GMPEs and GMICEs. In addition to this a model comparison will be made. ShakeMap models created by this thesis is compared with models by the USGS to highlight similarities and possible deviations of how the software is configured.
4.1.1 Scripts for running ShakeMap
One of the main objectives of this thesis is to create an earthquake intensity database that will be used by the IBEC classification scheme. This database is created by tabulating intensity data for a wide range of different earthquake scenarios modelled using ShakeMap. ShakeMap requires an input text file, with earthquake data, to create those models. For the purpose of modelling many scenarios at the same time, the process of editing these input files and converting the output to useful matrices may be tedious and time consuming. Therefore three separate Python scripts have been written that configure, run, and convert the ShakeMap required input files. The final output from these scripts are matrices used in the intensity approximation function. These scripts also simplifies future work as these scripts can be used to model other GMPEs than those modelled here, and to model scenarios with more site specific settings such as: Vs30, other magnitude and depth intervals, or different epicenter coordinates. Thus a possibility to go from global to local models. The individual purposes of the three scripts are:
• First : configures ShakeMap required input files with whichever magnitudes, depths and epicenter latitudes and longitudes and store these files in appropriate folders.
• Second : runs the input files configured in First through ShakeMap and gives all kinds of output that ShakeMap can produce. For example the intensity maps presented herein. Every map produced comes in text format as well, containing contour values for different coordinates.
• Third : converts the text files created in Second into intensity matrices which are used in the intensity approximation function.
The matrices given by Third are created by finding nearest intensity for a particular location in the vicinity of the earthquake epicenter. This process will be referred to as stepping and it is applied to the text file version of the intensity contour maps created by ShakeMap. The text files contain three columns where two of them are latitudes and longitudes, and the last one are corresponding intensity values. The stepping starts at the earthquake epicenter and extends to an epicentral radial distance of 555 km (∼5 degrees) and step size is kept constant at 4.66 km (∼0.042 degrees). For the global models the intensity contours are radially symmetric and the distance between epicenter and any contour is converted from latitude and longitude to kilometre to by-pass the fact that one unit of longitude is different depending on where you are on the globe. The stepping process is done for all modelled GMPEs and over a magnitude interval of 4.0 ≤ Mw≤ 9.0 with an increment of 0.5 Mwand a depth interval of 10 ≤ d ≤ 150 with an increment of 10 km. For each magnitude there will be 15 vectors, one for each hypocentral depth and these vectors are stored in a 2D matrix, one matrix for each magnitude. In total there will be 11 matrices of size [15 × 120] per GMPE.
4.1.2 GMPEs
As SNSN analyses tele-seismic earthquakes ≥ 5.0 Mw it is preferable to model GMPEs defined for earthquakes 4.0 ≤ Mw≤ 9.0. Though none of the available GMPEs are defined over exactly this interval.
What will be modelled and implemented into IBEC are GMPEs that are as close to this magnitude interval as possible and these are Akkar and Bommer (2007) (AkkarBoomer07), Akkar et al. (2013) (ASB13), Motazedian and Atkinson (2005) (MA2005) and Youngs et al. (1997) (Youngs97). See Table 1 for more information about each equation. However, only Youngs97, MA2005 and AkkarBoomer07 are discussed in detail and presented with graphs. How much the uncertainty of the intensity approximation increases while modelling outside defined intervals is not yet investigated. Note that all five GMPEs are available and operational in the intensity approximation function, mmiApprox(M, d, G, ¯r) (M: earthquake magnitude, d: source depth, G: chosen GMPE, ¯r: surficial distance to source). The user can simply change the value given to G to also change the GMPE used for approximating intensities.
The GMPE Youngs97 is chosen as the first choice when modelling ground motions. This means that the IBEC classification scheme and all scheme-related figures shown herein are based on this equation.
Youngs97 is a well-referenced GMPE (Wald et al. (2005)); it is designed using data from a subduction zone region and by analysing ShakeMaps created by the USGS, it is shown that they often use Youngs97 as well. Also, according to Ruff and Kanamori (2012) most high magnitude earthquakes come from subduction zones. Further, high magnitudes cause the greatest intensities and damages. Thereby, having a widely used subduction zone related GMPE, like Youngs97, when developing a global model for clas- sification of hazardous earthquakes is considered appropriate. How this assumption scale with normal
and transform fault earthquakes is yet to be concluded. The GMPE found by Youngs et al. (1997) is
ln(PGA)i j = C∗1+C2Mi+C3∗ln[(rrup)i j+ exp(C∗4C2
C∗3Mi)] +C5Zss+C8Zt+C9Hi+ ηi+ εi j
C∗1= C1+C3C4−C3∗C∗4
C3∗= C3+C6Zs
C4∗= C4+C7Zs
where i is the earthquake index, j is index of the station recording the ith event, M is Mw, rrup is distance between source and site in kilometres, H is source depth in kilometres, C are regression coef- ficients determined by the analysis and Z are indicator variables which takes the value 0 or 1 depending on the geological situation. For example, Table 3 summarizes regression coefficients of their spectral acceleration relationship. Spectral acceleration is used for, and a better estimator, studying motions of buildings. Although, normally spectral acceleration and peak ground acceleration do correlate. Figure 8 show how Youngs97 peak ground motions, converted to intensity, decrease in amplitude with increased distance to epicenter for various source depths. The figure reveals that intensity continuously decreases as source depth increases. For Mw≤ 8.0, the average intensity drop is 2.5 MMI over the whole depth in- terval making it somewhat possible to distinguish between shallow and deep earthquakes just by looking at intensity. The almost depth constant behaviour of intensity for Mw8.5 and 9.0 is interpreted represent the fact that earthquakes of such magnitudes rupture to 150 km depth sometimes. Or it could possibly be artefacts due to running models outside the GMPE’s intentional magnitude interval being 5.2 ≤ M ≤ 8.0.
The GMPE MA2005 is based on recordings done in the Puerto Rico region. Figure 9 show how MA2005 predict changes in peak ground motions, converted to intensity, and the structure of the figure is the same as for Youngs97. MA2005 seems to predict fairly constant values of intensities for sources starting at 40 to 70 km depth and downwards, for Mw≤ 7.0. Though intensity is a function of depth which is a known parameter and thus distinguishing deep from shallow earthquakes is still possible. Intensity is also dependent on seismic energy and the energy released by deep earthquakes is less concentrated at the surface than by shallow earthquakes due to geometrical spreading. Thereby, intensity is not likely to be constant with varying hypocentral depth. Though, earthquake sources are, as mentioned in Section 2.1, plane sources that span both space and time. A consequence of this is that the greatest movements along a plane is not necessarily in the hypocenter. A deeply located hypocenter may cause the same amplitude
Table 3. Some values of regression coefficients of an attenuation relationship for horizontal response spectral acceleration (5 % damping) for subduction earthquakes (modification from Youngs et al. (1997).
For rock
Period (s) C1 C2 C3 C∗4 C5∗
PGA 0.0 0.0 -2.552 1.45 -0.1
0.075 1.275 0.0 -2.707 1.45 -0.1
0.1 1.188 -0.0011 -2.655 1.45 -0.1
1.0 -1.736 -0.0064 -2.234 1.45 -0.1
1.5 -2.634 -0.0073 -2.160 1.50 -0.1
3.0 -4.511 -0.0089 -2.033 1.65 -0.1
For soil
Period (s) C1 C2 C3 C∗4 C5∗
PGA 0.0 0.0 -2.329 1.45 -0.1
0.075 2.400 -0.0019 -2.697 1.45 -0.1
0.1 2.516 -0.0019 -2.697 1.45 -0.1
1.0 -2.870 -0.0066 -1.785 1.45 -0.1
1.5 -5.101 -0.0114 -1.470 1.50 -0.1
3.0 -6.672 -0.0221 -1.347 1.65 -0.1
of ground motions as shallower sources which may explain the almost constant values of intensity seen in the figure.
Another considered GMPE is AkkarBommer07 by Akkar and Bommer (2007) and it is created using data recorded in Europe. Its intensity variations is seen in Figure 10. Figure 10e clearly reveals that AkkarBommer07 assumes intensity to be constant with depth. The author assumes that such behaviour is not optimal for earthquake classification. IBEC uses a global model, where incoming tele-seismic events can have any source depth and subduction zone generated earthquakes are often relatively deep.
For example, AkkarBommer07 assumes the intensity caused by a magnitude at great depth to be identical to a shallow earthquake of same magnitude. However, this GMPE could be, yet not investigated, useful for classifying earthquakes where source depth is unknown, less important or the earthquakes are always shallow.
(a) Source depth: 20 km (b) Source depth: 40 km
(c) Source depth: 60 km (d) Source depth: 80 km
(e) 0 km distance looking into the Earth
Figure 8. Variations of ground-motions predicted by Youngs97. Sub-figures (a) to (d) correspond to fixed source depths while (e) show variation of ground-motions at the epicenter and looking into the Earth.
(a) Source depth: 20 km (b) Source depth: 40 km
(c) Source depth: 60 km (d) Source depth: 80 km
(e) 0 km distance looking into the Earth
Figure 9. Variations of ground-motions predicted by MA2005. Sub-figures (a) to (d) correspond to fixed source depths while (e) show variation of ground-motions at the epicenter and looking into the Earth.
(a) Source depth: 20 km (b) Source depth: 40 km
(c) Source depth: 60 km (d) Source depth: 80 km
(e) 0 km distance looking into the Earth
Figure 10. Variations of ground-motions predicted by AkkarBommer07. Sub-figures (a) to (d) correspond to fixed source depths while (e) show variation of ground-motions at the epicenter and looking into the Earth.
4.1.3 GMICEs
Ground Motion/Intensity Conversion Equations (GMICEs) are used to convert peak ground motion (PGM) data to intensity values on the MMI scale (see previous Section 2.3). ShakeMap has two im- plemented GMICEs, namely Wald99 (Wald et al. (1999)) and WGRW11 (Worden et al. (2012)). This project uses Wald99, which was established specifically for ShakeMap (Wald et al. (2005)). The empir- ical relationship by Wald et al. (1999) for Modified Mercalli intensities 5 ≤ Imm≤ 8 for PGA is,
Imm= 3.66log(PGA) − 1.66 and for PGV for 5 ≤ Imm≤ 9,
Imm= 3.47log(PGV) + 2.35 and for Imm< 5,
Imm= 2.20log(PGA) + 1.00 Imm= 2.10log(PGV) + 3.40
They found their regression relationship by correlating peak ground motions to observed intensities for eight significant earthquakes occurring in California, U.S., between year 1971 to 1992 and earthquake magnitude 5.8 to 7.3 Mw. The correlating values of Imm are presented in Table 4 (refer to Table 2 for description of intensities).
Worden et al. (2012) found their relationship by correlating a database of approximately 200 000 in- tensity observations of California, U.S., earthquakes from USGS DYFI reports with peak ground motion amplitudes for earthquake magnitudes between 3.0 to 7.3 Mw. They claim their results to be in close agreement with the relationships of Wald et al. (1999) for Imm> 5. For lower intensities there is a slight deviation.
4.1.4 Comparison to USGS
Verification the ShakeMap installation and general configuration have been done by re-creating real events as scenarios and compare those with ShakeMaps created by the USGS for same real events. The author assumes that deviations of ±0.5 MMI should be considered fine but, more importantly, that the intensities modelled herein using ShakeMap seems to follow the topography in the same way as the USGS models.
Table 4. Ranges of peak ground motions for Modified Mercalli Intensities (from Wald et al. (1999)).
Intensity Peak Acceleration
(% g) Peak Velocity (cm/s)
I < 0.17 < 0.1
II-III 0.17 - 1.4 0.1 - 1.1
IV 1.4 - 3.9 1.1 - 3.4
V 3.9 - 9.2 3.4 - 8.1
VI 9.2 - 18 8.1 - 16
VII 18 - 34 16 - 31
VIII 34 - 65 31 - 60
IX 65 - 124 60 - 116
X+ > 124 > 116
Figure 11 and Figure 12 show two different events: Mw7.1 at 40 km depth that occurred in Peru and Mw6.1 at 28 km depth that occurred in New Zealand. Each figure contains two sub-figures, where the left sub-figure is created by this project and the right sub-figure is created by the USGS. Both the Peru and New Zealand earthquakes modelled by this project uses Youngs97 as their GMPE. This will give an insight into how the result of using only that equation scales with the results achieved by other GMPEs.
The USGS models the Peru earthquake using Youngs97 and the New Zealand earthquake using Zhao06.
Comparing Figure 12a with 12b one sees that the ShakeMap configuration used within this project sug- gests a slightly higher intensity than what is suggested by the USGS while Figure 11 shows the opposite and this configuration suggests a slightly lower intensity. These differences are thought to appear as this project uses the GMICE Wald99, while the USGS uses the GMICE WGRW11, for both events. An- other reason could be differences in extracting the Vs30data and minor configuration differences within ShakeMap. For this comparison we used a Vs30data generator (USGS (2015a)) provided by the USGS and assumes that Peru is an active tectonic region rather than a stable shield. How the USGS themselves generate their Vs30is unknown.
(a) ShakeMap by IBEC (b) ShakeMap by USGS
Figure 11. Verifying the ShakeMap installation by re-creating real earthquakes and compare the same earthquake model done by the USGS. Mw7.1 at a hypocentral depth of 40 km.
(a) ShakeMap by IBEC (b) ShakeMap by USGS
Figure 12. Verifying the ShakeMap installation by re-creating real earthquakes and compare the same earthquake model done by the USGS. Mw6.1 at a hypocentral depth of 28 km.
4.2 The SNSN classification scheme
For earthquake classification the SNSN scheme uses three classes: A = high damage, B = medium damage and C = low potential damage. To decide which class to assign the scheme uses the following parameters:
• earthquake magnitude
• hypocentral depth
• epicenter latitude and longitude
• population density
• distance to coast line
• elevation (negative for bathymetry)
• in ocean or on land
• closest plate boundary
The importance of earthquake magnitude and hypocentral depth is trivial and greatly controls how intense the ground motions are. With the use of epicenter latitude and longitude the SNSN scheme can extract additional data related to geography and the society such as population density and closest plate boundary. The latter is the deciding factor of whether a possible tsunami is generated or not, as tsunamis require vertical movements along a plate boundary. When estimating the size of population density affected by the earthquake the SNSN uses a search radius, over which the density is calculated, that is equivalent to an average radial fault length as a function of magnitude. For instance, Mw≥ 8 estimates R to be 55 km, while Mw< 6 estimates R to be 5 km.
4.3 The Intensity Based Earthquake Classification (IBEC) scheme
Regarding the run-time for IBEC (Figure 13) there is no apparent difference compared with the SNSN scheme, even if the IBEC scheme provides more information. The reason is because all intensity data are pre-defined and stored in separate, relatively small, matrices and once an incoming earthquake is about to be classified it does not require the whole intensity database to be opened. A complete run for one event is about 2 seconds where most time is spent on extracting and estimating the population density.
The IBEC method estimates the search radius, used for estimating population density, as the radial distance to an intensity contour. Used distance have a minimum and maximum of 15 and 555 km, respectively. Which contour depends on the epicentral earthquake intensity, I0. If I0< 7.0, search radius is the radial distance of how far away an intensity of at least 5.5 MMI will affect the ground. If I0≥ 7.0, the function instead returns the distance to the 6.5 MMI contour. The argumentation behind this change is: if the predicted epicentral intensity is high, severe earthquake damage is more likely to occur and thus easier to predict. The methodology and its results is further discussed in subsequent sections.
Magnitude, depth and epicenter coordinates
Intensity
GDP per capita MMI
Pop. dens.
Pop. dens.
C B A
MMI
Pop. dens. Pop. dens.
C B
A
mmiApprox()
≤20000
> 4.0
≤ 4.0
≤ 100
any > 100
>
20000
> 6.0 ≤ 6.0
≤ 300 any
> 300
Figure 13. Flowchart of a generalized version of the IBEC classification scheme. The boundary values are not necessarily the values used in the real version, but they are not out of context and should give an intuitive feeling of how the classification works. A, B and C are the classes assigned to an earthquake and assumes the following of death tolls: A ≥ 100, 1 ≤ B < 100, and C = 0.
4.3.1 Earthquake intensity
The intensity matrices explained in Section 4.1.1 have enabled the possibility to define the function mmiApprox() (see Appendix B for the function’s code). It gives the earthquake intensity experienced at locations in vicinity of the epicenter. Such locations could be nearby cities, nuclear power plants and other important facilities. Intensity varies between 1 to 10 MMI. The function requires earthquake magnitude and source depth to be known and also a GMPE from Table 1 has to be chosen. The algorithm flowchart is:
Magnitude, M, and depth, d
Locations, ¯r, and GMPE, G
Load intensity matrix, A(M, G)
Intensity = A[row(d), col(¯r)]
mmiApprox(M, d, G, ¯r) is independent of epicenter’s latitude and longitude as the intensity matrices are modelled without any site-specific settings, such as topography and geology. The vector ¯r is created using the Haversine formula, haversineFormula() (see Appendix D for the function’s code), which cal- culates the distance in kilometre between two locations on a sphere, given their latitudes and longitudes which are in DD.MMMM format. The Haversine formula assumes that the Earth is perfectly spherical with an average radius of 6371 km. The real radius at the poles and equator are 6357 km and 6378 km, respectively.
All events from the SNSN database have been ran through the mmiApprox() function and plotted in a magnitude-depth-epicentral intensity environment. The result is seen in Figure 14. Figure 14a plots all events being unfiltered while Figure 14b plots all events that have occurred in areas with a population density > 0. The area enclosed by the red dashed lines depicts the magnitude and depth intervals used for modelling scenarios to create the intensity matrices and is the true area over which mmiApprox() is defined, 4.0 ≤ M ≤ 9.0, 10 ≤ d ≤ 150. If the function is ran with magnitudes or depths outside this area these parameters will be rounded off to the closest value within the area. For example, an earthquake of Mw3.4 at 8 km depth will be treated as a Mw4.0 at 10 km depth while a magnitude Mw 9.2 at 200 km depth will be treated as a Mw9.0 at 150 km depth. Note the group of markers at 550 to 650 km depth in Figure 14b that appears when filtering on events with a population density > 0. This pattern is thought to appear as subducted plates tend to rest at the 670-discontinuity and deeply subducted plates have dived in under the continents, which then are populated. This leaves us with a group of earthquakes with depths of 300 to 550 km occurring under areas of zero population density. These earthquakes could be related to non-populated volcanic arcs such as the Andes, and back arc regions such as the ocean between Japan and China.
(a) All events. (b) Events with population density > 0.
Figure 14. Earthquakes from the SNSN data base have been ran through mmiApprox() using the GMPE Youngs97.
The red dashed lines represent the magnitude and depth intervals used for modelling scenarios to create the intens- ity matrices. Earthquake intensity increases with magnitude, and decreases with depth.
The circular surface area over which population density is estimated is related to the area which is affected by ground motions. The SNSN scheme uses a radius that is a function of earthquake magnitude.
The IBEC suggests and uses a method where radius is a function of earthquake magnitude, depth and intensity. The IBEC method could be explained by imagining the Earth’s surface above a hypocentrum being a grid where each grid node is an earthquake intensity. The size of the grid is determined by a minimum intensity. A suitable minimum would be the lowest intensity that have been reported to cause casualties. The IBEC uses a minimum of 5.5 MMI and the corresponding grid for two different events are the yellow areas in Figure 15. It is seen that different earthquakes can cause the same intensity at the surface but over completely different areas. Figure 16 shows the distance on the surface to the MMI 5.5 contour at different depths for the Mw7.5 earthquake seen in the previous figure. It is seen that distance increases for sources at 10 to 110 km depth. Beyond 110 km depth the distance instead decreases because of the progressively stronger attenuation and the energy density and concentration decreases as source depth increases. Thereby, by considering magnitude, depth and intensity altogether we are able to make a good estimate of the size of the area that will shake remarkably. The function written that makes these estimations is distToCont() (Appendix C) and it assumes a flat surface, and does not take topography into account. The function’s flowchart is:
Magnitude, M, and depth, d
Contour, c, and GMPE, G
Load matrix, A(M, G)
Index of A[row(d), col(c)]
Radius
= D[index]
The algorithm requires four parameters: magnitude, M, source depth, d, GMPE, G and intensity, c, which is the intensity that the function will calculate the epicentral distance to. The magnitude together with the chosen GMPE defines a matrix where each matrix entry represents an intensity value. The algorithm then finds the matrix entry that is equal to c. When the decrease in earthquake intensity with distance is slow there are more than one matrix entry that equals c. In those cases the function picks the entry that is furthest away from the epicenter. Once the correct entry is found it takes its index. This index is then put into a separate distance vector, D, that is of same length as the matrix rows. Thus putting the entry index into D returns the distance in kilometre from the epicenter to the position on the surface that is experiencing an earthquake intensity of c. In extreme cases, however, where the last entry of the intensity row given by A[row(d), :] is larger than c, c is updated to instead be the value last entry of A[row(d), :]. From this follows that the search radius will be the very last entry of the position vector that is 555 km. Such extreme cases arise, for example, if the algorithm is fed with with magnitudes around 9.0 and depths around 50 km while looking for an intensity of 5.5 MMI or lower.
(a) Mw6.5 at 10 km depth. (b) Mw7.5 at 110 km depth.
Figure 15. Variations in size of surface area covered by an earthquake intensity of, in this case 6.0 (yellow color), for two earthquakes of different magnitudes and depths.
Figure 16. Graph of how the surface distance from epicenter to the intensity contour of MMI 5.5 varies with source depth for a Mw7.5 earthquake. The GMPE used is Youngs97.
4.3.2 Near-shore earthquakes
Figure 17 is an illustration of how near-shore earthquakes are treated in the classification process. R is the radial distance to an user-defined intensity, and it is now set to 5.5 MMI. dcis the distance to nearest coast line. If R − dc< 0 the defined intensity never reaches the coast and the land will not experience ground motions high enough for the earthquake to get a classification other than C. Instead if R − dc ≥ 0 the defined intensity overlap with the nearest coast line and the classification will cover all three classes, A, B and C, and be determined by the size of the population in the area and earthquake intensity. A shortened version of the generalized classification scheme for near-shore earthquakes is illustrated by Figure 18.
Its flow chart looks like the generalized IBEC scheme in Figure 13 apart from that a ’if R - dc <
0, elif R - dc >= 0’ statement would take place right before the GDP per capita determination.
The continuation of the scheme after the node ’GDP per capita’ is the same as in Figure 13, but is not shown in this flow chart to prevent unnecessary repetition.
Figure 17. Illustration of the parameters introduced in IBEC when classifying a near-shore earthquake: epicenter’s distance to nearest coast line, dc, and distance from the epicenter to a pre-defined earthquake intensity, R, also referred to as search radius. The star is the earthquake epicenter.
Intensity
Contour- versus coast-distance
GDP per capita C
R−dc
<0
R− dc≥ 0
Figure 18. Flowchart of how near-shore events are treated by IBEC in the node Contour- versus coast-distance.
Note that this is a downsized version and only the necessary parts, regarding near-shore earthquakes, of the classi- fication scheme are shown. The C class correspond to zero death tolls.
4.3.3 Building standard
Building infrastructure and houses that are resistant to ground motions, preventing them from a collapse, is a good way of reducing the risk of an earthquake causing casualties. How these buildings should be built is specified by so-called building codes which are established by engineers studying, for example, simple harmonic oscillators. Such studies require funding and the buildings themselves often become more expensive to construct if codes are used. With that said, an earthquake with an intensity of 7.0 MMI in one location might not cause any damage, while the same intensity somewhere else might cause a total disaster. Knowledge of whether building codes are used or not is useful when a scheme is classifying
earthquakes and predicting their damages. However, databases of if, or to what degree, countries are using building codes do not exist or are not consistent. On the contrary, databases of GDP per capita are easily found and globally defined. The IBEC scheme assume GDP per capita to be an indicator of whether building codes are used or not. A high GDP per capita suggests a high quality of buildings. A low GDP per capita suggests the opposite, a low quality of buildings.
4.3.4 Population
The Modified Mercalli Intensity scale says that at intensities of 5.0 to 6.0 MMI windows and plaster starts to break and thus becomes a threat to humans. Also, by analysing the death toll data described in Section 3.1 it is seen that death toll increases from zero at intensities of 5.5 MMI. This value of 5.5 is used in the function for determining which search radius, R, that is used when estimating earthquake affected population densities.
An intensity of MMI 5.5 or above also works as a threshold for when IBEC evaluates if estimated population density is not reasonable high relative to search radius. One could argue that the larger the search radius the higher the population density. But if an epicenter is in a remote area such as deserts and mountainous areas with only one city the population density for the whole affected area will be low. That is, the second closest city lies outside the MMI 5.5 contour. So when these thresholds are reached IBEC classifies the earthquake using an integrated population of individual affected cities and the intensity experienced in those cities, instead of epicenter intensity. For example, a Mw 6.0 at some depth will have a search radius of around 60 km and the area affected by an intensity of 5.5 MMI is then 2πR2≈ 23000 km2. Further assume that the earthquake occurs in a remote area with only one city, and neighbouring cities are located far away. If this earthquake occurs right beneath the populated city then the population density will be low compared to the large search area given above. The city is at high risk of getting damaged but a general classification scheme will not realise this, and will not classify this earthquake as a high damage potential earthquake because of low population density. But IBEC compares the intensity experienced in cities surrounding the epicenter with the overall population density to better classify earthquakes occurring in remote areas.
4.3.5 Case studies
In this section we explain how the intensity approximation function works using an arbitrary earth- quake. Imagine an earthquake of Mw 7.4 that occurs at a hypocentral depth of 34 km, with coordinates 5◦15’56.3” N and 95◦37’46.8” E, and the GMPE Youngs97 has been pre-defined by the user. Using this information, intensity at epicenter is Imm= 7.0 MMI. The SNSN function FindCities(lat, lon, ...) reveals that the earthquake occurred in Indonesia and thus GDP per capita is found to be 10200 Int$. Thereby,
this particular event will take the left route out from the node ’GDP per capita’ in Figure 13. Further on, as Imm= 7.0 was determined, the event will take the path straight down as it fulfils the requirement Imm> 4.0. Arriving at the new node, ’Pop. dens.’, the population density affected by an intensity ≥ 5.5 will be determined by first calculating the distance from the epicenter to the 5.5 intensity contour using R= 154 km. Population density is found by using another SNSN function p = 2614 humans/km2. Now with all necessary data given, the classification of this particular event is A which is the highest class and it predicts a death toll ≥ 100. The other classes are B and C which predict death tolls of 1 → 100 and 0.
There are also two big nearby cities named Banda Aceh and Lhokseumawe with coordinates 5◦33’0”
N, 95◦19’0” E and 5◦11’17” N, 97◦8’25” E, respectively. Due to their importance, it is necessary to more specifically approximate the intensities in these cities. This is done by calculating the distance from the epicenter to each city using ¯r = [46, 167] km. The resulting earthquake intensities in Banda Aceh and Lhokseumawe are 7.0 and 5.0 MMI, respectively. At this point IBEC does not make separate classifications of the potential damage in each individual city unless they are in epicentrum. The scheme classifies the overall potential damage of an earthquake.
