i SEISMIC HAZARD ASSESSMENT AND VOLCANOGENIC SEISMICITY FOR THE DEMOCRATIC REPUBLIC OF CONGO AND SURROUNDING AREAS, WESTERN RIFT VALLEY OF AFRICA Georges Mavonga Tuluka A thesis submitted to the Faculty of Science, University of the Witwatersrand, in fulfilment of the requirements for the degree of Doctor of Philosophy Johannesburg, 2010 ii DECLARATION I declare that this thesis is my own, unaided work. It is being submitted for the Degree of Doctor of Philosophy in the University of the Witwatersrand, Johannesburg. It has not been submitted before for any degree or examination in any other University Georges Mavonga Tuluka On the day of 20 iii ABSTRACT The Western Rift Valley of Africa has experienced several severe earthquakes with magnitude M?6 and volcanic eruptions in recent historical times. Most of these earthquakes occurred in the Democratic Republic of Congo (DRC) and neighbouring countries such as Uganda and Tanzania. The largest earthquake on record, which occurred at Kasanga (Tanzania) on 13 December 1910 in the Lake Tanganyika area, had a magnitude of Ms7.3. The most recent significant earthquake occurred approximately 20 km north of Bukavu City (DRC) on 03 February 2008 with magnitude Mw=6.0 in the Basin of Lake Kivu. The Virunga volcano group, located at the northern edge of Lake Kivu, consists of eight major volcanoes (Muhavura, Gahinga, Sabinyo, Visoke, Karisimbi, Mikeno, Nyiragongo and Nyamuragira). The volcanoes Nyiragongo and Nyamuragira have been the most active since 1882. A probabilistic approach was used to map the seismic hazard in Democratic Republic of Congo and surrounding areas, and assess the seismic hazard level for 14 cities in the region. Seismic hazard maps for 2%, 5% and 10% chance of exceedance of the indicated ground accelerations in 50 years were prepared using a 90-year catalogue compiled for homogeneous magnitudes (Mw); the attenuation relations of Mavonga (for the Western Rift Valley of Africa), Atkinson and Boore (for eastern and North America) and Jonathan (for eastern and southern Africa); and the EZ-Frisk software package. A Poisson model of earthquake occurrence that assumes that events are independent was adopted. Therefore, foreshocks, aftershocks and earthquake swarms were removed from the initial catalogue of 2249 events. Furthermore Mw=4 was selected as the lower magnitude bound (Mmin) because smaller earthquakes are considered unlikely to cause damage, even to houses that are poorly designed and built. Thus, any remaining events with Mw<4 were also excluded from the catalogue, leaving a sub-catalogue of 821 events iv The highest estimated levels of seismic hazard were found in the Lake Tanganyika Rift seismic zone, where peak ground accelerations (PGA) in excess of 0.32g, 0.22g and 0.16g are expected to occur with 2%, 5% and 10% chance of exceedence in 50 years, respectively. The seismic hazard in the Congo basin diminishes with distance away from the Western Rift Valley until, at a distance of about 450 km, the chance of exceeding 0.05 g (the threshold value of engineering interest) is less than 10% in 50 years. Finally, from the probabilistic seismic hazard analysis of the DR Congo and surrounding areas, four seismic zones were identified and rated as follows: Zone A (very high hazard), Zone B (high hazard), Zone C ( moderate hazard), and Zone D (low hazard). The zone A includes the Lake Tanganyika Rift zone where PGA values of 0.32g, 0.22g and 0.16 g are expected to occur with probability 2%, 5% and 10% in 50 years, respectively. Zone B includes the Lake Kivu Basin, Mount Ruwenzori and Lake Edouard region. Zone C includes Rutsuru, Masisi, Upemba area and a part of the Congo basin close to the Western Rift. The remainder of the Congo basin constitutes the zone D. To understand how volcanoes work and reduce the risk due to the Virunga volcanic eruptions in the Western Rift Valley of Africa, the crustal structure beneath the Virunga volcanic area was investigated and studies of volcanogenic seismicity were carried out. From these studies, it was found: High velocity material (6.9 to 7.3 km/s) lies beneath the Kunene (KNN) and Kibumba (KBB) broadband stations at depths of 3-20 km and 3- 10 km, respectively, which is indicative of magma cumulates within the volcanic edifice. A low velocity zone was observed below KNN and KBB at depths of 20- 30 km and 18-28 km, respectively, and with average velocity 6.1 km/s and v 5.9 km/s. This low velocity zone may sample the magma chamber or a melt-rich sill. The depth of the crust-mantle transition zone beneath the area sampled by the KNN and KBB in the Virunga area was determined to be about 39 to 43 km and 30 to 39 km, respectively. Swarm-type seismicity composed mainly of long-period volcanic earthquakes preceded the eruptions of Nyamuragira volcano and was probably enhanced by tectonic seismicity related to rifting. A steady increase in seismicity at a constant rate from a deep magma feeder (located at about 20 to 30 km depth) was observed ten or eleven months before eruption. In the last stage (1 or 2 months) before the eruption, the hypocenters of long-period volcanic earthquakes became shallower. The new model of the local crustal seismic velocity for the Virunga area is useful to map the migration of hypocenters of earthquakes accurately and reveals trends that could be precursors of volcanic eruptions This pattern of seismicity prior to the volcanic eruptions, integrated with other available data (e.g. INSAR, GPS), may be used to characterize the volcanic process and forecast volcanic eruptions in the Virunga area. vi ACKNOWLEDGEMENTS This study was carried out at University of the Witwatersrand, South Africa with full financial support from the Belgian Technical Cooperation. I am sincerely grateful to Belgian Government. A number of individuals have provided support in different aspects and stimulated me in my work. I owe my sincere thanks to Professor Ray Durrheim (my supervisor) for his wise guidance and encouragement during my work, to Professor A. Kijko (Aon Benfield Hazard Centre, Univeristy of Pretoria, South Africa), for his suggestions and documentations during the initial stages of this study. I am indebted to Emeritus Professor Hiroyuki Hamaguchi and Associate Professor Takeshi Nishimura of Tohuku University, Japan, for installing digital seismographs (SEIDAS) in the Kivu Province, DR Congo which allowed me to get quality data to study major earthquakes which occurred from 1994 to 2002 and determine attenuation relationship for strong motion in the Kivu Province. To Professor Zana Ndontoni, our former General Director of ?Centre de Recherche en Sciences Naturelles de Lwiro? and acting General Director of ?Centre de Recherche en Geophysique de Kinshasa?, for his encouragement and giving young people opportunities to study seismology. I owe him many thanks. My sincere appreciations are due to co-Directors of AfricaArray, Professors Andy Nyblade and Paul Dirks and others for their initiative of launching this program which is primary educational initiative to strengthen geosciences training and research program in support of Africa?s natural resource sector and geohazard vii assessment in Africa. To Dalena Blitenthall, Melody van Wyngaard and all staff of the Secretariat of AfricaArray and the School of Geosciences at Witwatersrand University, my sincere thanks for assisting in many administrative tasks and logistics. Your warm friendly attitude and assistance created a stimulating atmosphere that was conductive to serious work. I further owe my thanks to all my colleagues and staff of Goma Volcano Observatory for assistance with data collection, Management Committee of ?Centre de Recherche en Sciences Naturelles de Lwiro? and the Ministry of Education and Scientific Research of DR Congo, for granting me a part of my salary at times I was away from my office. I am sincerely grateful. Many thanks to my wife Annie Lukombo and my children (Divine, Lydia, Chrislin, Dieu Merci, Guershom and Obed), for their patience, love, encouragement and moral support. Last but not the least; I thank all those, even though not mentioned by name that stimulated me to study seismology. viii CONTENTS DECLARATION......................................................................................................i ABSTRACT............................................................................................................iii ACKNOWLEDGEMENTS....................................................................................vi LIST OF FIGURES................................................................................................xi LIST OF TABLES ................................................................................................. xv GLOSSARY OF TERMS .................................................................................... xvi LIST OF OWN PAPERS.....................................................................................xxv CHAPTER 1 INTRODUCTION 1.1 General Introduction .......................................................................................... 1 1.2 Motivation for this study .................................................................................... 5 1.3. Aim and structure of the study .......................................................................... 7 CHAPTER 2 GEOLOGY AND SEISMOTECTONICS OF THE DRC AND ADJACENT AREAS 2.1 Introduction ...................................................................................................... 11 2.2 Main seismic source zones for the study area .................................................. 12 2.2.1 The Western Rift Valley of Africa (WRA) .......................................... 12 2.2.1.1 Southern Sudan, Ruwenzori area, and Lake Edouard trough ... 14 2.2.1.2 Virunga volcanic complex, Rutsuru and Masisi area ............... 15 2.2.1.3 Lake Kivu Basin, Ngweshe area and Ruzizi plain ................... 16 2.2.1.4 Lake Tanganyika Rift ............................................................... 19 2.2.2 Southeastern DRC and northwestern Zambia ....................................... 19 2.2.3 The Congo basin ................................................................................... 24 ix CHAPTER 3 EARTHQUAKE HAZARD ASSESSMENT FOR DRC AND SURROUNDING AREAS 3.1 Data .................................................................................................................. 26 3.2 Method ........................................................................................................... 288 3.2.1 Unification of magnitudes .................................................................. 288 3.2.2 Completeness of the catalogue ............................................................. 32 3.2.3 Probabilistic seismic hazard analysis .................................................... 34 3.2.4 Seismic source zones ............................................................................ 35 3.2.5 Seismic parameters of seismic source zones ........................................ 37 3.2.5.1 b-value and activity rate ? ........................................................ 38 3.2.5.2 Maximum magnitudes .............................................................. 40 3.2.5.3 Focal depths .............................................................................. 43 3.2.5.4 Attenuation relations ................................................................ 43 3.2.6 Input parameters for hazard computation ............................................. 49 3.2.7 Computation of seismic hazard ............................................................ 50 3.2.7.1 Theoretical background ............................................................ 50 3.2.7.2 Seismic hazard analysis ............................................................ 51 3.3 Results .............................................................................................................. 53 3.3.1 Seismic hazard map .............................................................................. 53 3.3.2 Total hazard curve at the sites of the cities ........................................... 57 3.4 Discussion ...................................................................................................... 67 3.4.1 General overview .................................................................................. 67 3.4.2 Source of uncertainties ......................................................................... 68 3.4.3 Assessment of the sensitivity of the outputs to input data .................... 69 3.4.4 Comparison with previous work in the Western Rift Valley ............... 71 CHAPTER 4 VOLCANOGENIC SEISMICITY IN THE VIRUNGA AREA 4.1 Introduction ...................................................................................................... 73 4.2 Data acquisition ................................................................................................ 77 x 4.3 Crustal structure beneath two seismic broadband stations in the Virungavolcanic area ............................................................................................. 78 4.3.1 Method .................................................................................................. 78 4.3.1.1 Receiver function ...................................................................... 78 4.3.1.2. Data processing ....................................................................... 82 4.3.2. Results .................................................................................................. 91 4.3.2.1 Ps converted phases and Moho depth ....................................... 91 4.3.2.2. Inversion results of average receiver function at KNN and KBB .............................................................................................................. 92 4.3.3. Discussion ............................................................................................ 97 4.3.3.1 Lateral variation in receiver function with azimuth ................. 97 4.3.3.2 Crust-mantle boundary beneath KNN and KBB .................... 100 4.3.3.3 High and low velocity zone .................................................... 101 4.3.3.4 Application of the new velocity model to the hypocenter determination ...................................................................................... 102 4.4 Seismic activity prior to the recent eruptions of volcano Nyamuragira ......... 104 4.4.1 Method .............................................................................................. 104 4.4.2. Temporal and spatial variations in seismicity ................................... 108 4. 4.3 Summary of observations .................................................................. 109 CHAPTER 5 CONCLUSIONS AND RECOMMENDATIONS 5.1. Conclusions ................................................................................................... 114 5.2. Recommendations ......................................................................................... 116 REFERENCES ..................................................................................................... 118 APPENDIX A ...................................................................................................... 137 xi LIST OF FIGURES Figure1.1 Large earthquakes with magnitude M?6 observed in the Western Rift Valley of Africa since 1910. ..................................................................................3 Figure 1.2 Damage to a house resulting from earthquake on 3/2/2008 at Bukavu City (D.R Congo). Source: I. Badriyo, Project UNOPS ?Unite de gestion de risque volcaniques?, Goma,North- Kivu, Democratic Republic of Congo..........?.........4 Figure 2.1 Hypsographic DEM of the East African Rift system. Black lines are faults. Source: Chorowicz, J., 2005....................................................................?13 Figure 2.2 Generalised geology map of the Kivu Province. Source: Modified from Twesigomwe (1997)........????????????.14 Figure 2.3 Map showing the eight volcanoes of the Virunga Volcanic field and seismicity of Lake Kivu and surrounding area for the period 1976-2002. [Modified from Komorowsky et al., 2003 and Mavonga, 2007a]...?????..17 Figure 2.4 Epicenter map for the period 1955 to 1962 showing the central part of the Western Rift Valley of Africa ( Masisi, Lake Kivu, Ngweshe and Ruzizi plain). Source: J. Wohlenberg, 1968].....................................................................18 Figure 2.5 A map of Eastern African Rift System (EARS) showing the continental breakup in the south East of DR Congo and northern Zambia [Modified from Kinabo et al., 2007].................................................................?.22 Figure 2.6 Regional gravity anomaly and heat flow determination values at the sites reported from Zambia by Chapman and Pollack [1977], Kenya and Tanzania by Nyblade et al. [1990] and from D.R. Congo by Sebagenzi et al. [1993]. [Modified from Sebagenzi and Kaputo, 2002]?????????????..23 Figure 2.7 (a) Map showing the regional setting of the Congo Basin (b) Focal mechanism plots of all studied earthquakes for source dynamics in the Congo Basin. [After Atalay, 2002]...............................................................?????.25 Figure 3.1 Regression of Ms versus Mb ISC?s magnitudes for events observed in the Western Rift Valley of Africa......................................................................29 Figure 3.2 Regression of M(LWI) versus Mb(USGS).........................................30 xii Figure 3.3 Frequency distribution of M?2 earthquakes in the 1910-2007 catalogue. Magnitude increment is (a) ?m=0.5 and (b) ?m=0.25. .......................33 Figure 3.4 (S1) Epicentres in the D. R. Congo and surrounding areas for the time period 1910-2007 used in the seismic hazard analysis. (S2) The seismic source areas are marked as polygons ...............................................................................36 Figure 3.5 Cumulative magnitude-frequency plot for a magnitude increment ?m= 0.25 for the three main seismic source zones used in seismic hazard analysis: (a) Upemba-Moero Rift, (b) Congo basin, and (c) Western Rift Valley..???.39 Figure 3.6 Attenuation relations curves obtained for Central and Southern Africa data, together with the relations for Eastern North America (Atkinson and Boore, 2006) and Central USA (Somerville et al., 2001).. ...............................................48 Figure 3.7 (a) Distribution of mean PGA values (in unit g) in the DR Congo and surrounding areas computed for 2 % chance of exceedance in 50 years. ............54 Figure 3.7 (b) Distribution of mean PGA values (in unit g) in the DR Congo and surrounding areas computed for 5 % chance of exceedance in 50 years.???..55 Figure 3.7 (c) Distribution of mean PGA values (in unit g) in the DR Congo and surrounding areas computed for 10 % chance of exceedance in 50 years...??..56 Figure 3.8 Total hazard curve at the site (a) Bujumbura (b) Bukavu...............?59 Figure 3.8 Total hazard curve at the site (c) Bunia (d) Butembo.........................60 Figure 3.8 Total hazard curve at the site (e) Goma (f) Kalemi.............................61 Figure 3.8 Total hazard curve at the site (g) Kananga (h) Kigali...?????..62 Figure 3.8 Total hazard curve at the site (i) Kigoma (j) Kindu........................?63 Figure 3.8 Total hazard curve at the site (k) Kisangani (l) Lubumbashi..???64 Figure 3.8 Total hazard curve at the site (m) Mbuji Mayi (n) Uvira.....???...65 Figure 3.9 Contribution of hazard by seismic source at Bukavu site....................66 Figure 4.1 Map showing the eight volcanoes of the Virunga volcanic zone....... 74 xiii Figure 4.2 Epicentre distributions of earthquakes used to determine P receiver functions.........................................................?????????????....83 Figure 4.3 Stacked receiver functions of waveforms obtained at (A) KNN and (B) KBB stations in the NE quadrant at azimuthal interval between 0? to 45?, 60? to 90? (Left) and 0? to 90? (Right)...????????????????.....87 Figure 4.4 Comparison of receiver functions of deep earthquakes and shallow events at azimuthal interval 90? to 105? in the SE quadrant at (A) KNN and (B) KBB stations...................................................................................???.....88 Figure 4.5 (A) Stacked functions representative of NE, SE, SW and NW quadrants along with the average receiver function (AVG) for all quadrants computed at KNN station....................................................................???.....89 Figure 4.5 (B) Stacked functions representative of NE, SE, SW and NW quadrants along with the average receiver function (AVG) for all quadrants computed at KBB station...................................................................???.......90 Figure 4.6 Synthetic waveforms produced by the three initial models which compare well to the observed average receiver function at (A) KNN and (B) KBB.................................................................................................................94 Figure 4.7 (Left) Acceptable velocity models obtained from the three initial models and (Right) average velocity model at (A) KNN and (B) KBB ??.......95 Figure 4.8 (A) Estimate of the crustal structure for the Virunga area. (B) Simplified crustal structure in the Virunga area..............................................96 Figure 4.8C A satellite image of the Virunga volcano group and Lake Kivu. The filed triangles and circles indicate seismographic stations and Moho piercing points of rays at 36 km depth of teleseismic events used for receiver function analysis................................................................................................................98 Figure 4.9: (a) Map of Bouguer gravity anomalies (5mgal) contour interval in the Virunga area (b) Two dimensional model of Bouguer anomalies observed along D-D? shown above density model and topographic profile. Source: Zana et al. [1992b].........................................................??????...99 Figure 4.10 Hypocenters map of volcanic earthquakes for the period from 26 November 2006 to 27 November 2006 computed using (A) GVO velocity model and (B) The new velocity model...............................................................103 xiv Figure 4.11 Velocity spectra of typical LP event (18 March 2004) and hybrid (LP+short-period) event (27 November 2006)....................................................106 Figure 4.12 Velocity models used in this study.................................................107 Figure.4.13 Monthly count of locatable (series 1) and located (series 2) long- period volcanic earthquakes 11 months before the 8 May 2004 eruption of Nyamuragira (Top) and the 27 November 2006 eruption of Nyamuragira (Bottom)...............................................................................................................110 Figure 4.14 Long and short-period earthquake hypocenters observed in the Virunga area during the period (A) 24 November 2002 to 24 December 2002 and (B) 01 July 2004 to 01 August 2004.............................................................111 Figure 4.15 Long and short-period earthquake hypocentres observed in the Virunga area during the period (A) 07 January 2004 to 07 February 2004 and (B) 25 July 2006 to 25 August 2006....................................................................112 Figure 4.16 Long and short-period earthquake hypocentres observed in the Virunga area during the period (A) 07 April 2004 to 07 May 2004 and (B) 28 October 2006 to 27 November 2006.........................................................113 xv LIST OF TABLES Table 3.1 Seismicity parameters for the area source zones...................................49 Table 3.2 Peak ground acceleration at selected cities determined using probabilistic and deterministic approaches for two seismic source zone models of the Western Rift Valley.....................................................................................58 Table 4.1 Eight major volcanoes of the Virunga area...........................................75 Table 4.2 Ps-P delay time and Moho depth..........................................................91 Table 4.3 Average local velocity model in the Virunga area................................93 Table A.1 List of earthquakes used in this study for the regression analysis of Mb (USGS) versus M (LWI)...........................................................................137 Table A.2 List of hypocenters and magnitudes of earthquakes used in the determination of receiver functions.....................................................................140 xvi GLOSSARY OF TERMS Source: Keiti, A. and Lee, W.H.K., 2003. Glossary of Interest to Earthquake and Engineering Seismologists. In Lee, H.K., Kanamori, H., Paul, C.J., Kisslinger, C. (editors), International Handbook of Earthquake & Engineering Seismology, Academic Press. Part B, Vol. 81 B, p. 1793- 1856, Academic Press, San Diego, California, USA Acceleration: The rate of change of particle velocity per unit time. Active fault: A fault that has moved in historic (e.g., past 10 000 years) or recent geological time (e.g., past 500 000 years). Although faults that move in earthquakes today are active, not all active faults generate earthquakes as some faults are capable of moving aseismically. Active volcano: By a widely used but poor definition, one that is erupting or that has erupted one or more times in recorded history Aftershock: An earthquake occurring as a consequence of a larger earthquake in roughly the same location. The sequence of such earthquake following a larger one generally shows a regular decrease in the rate of occurrence, first discovered by Omori in 1894. Aleatory uncertainty: The uncertainty in seismic hazard analysis due to inherent random variability of the quantity being measured. Aleatoric uncertain- ties cannot be reduced by refining modelling or analytical techniques. See also epistemic uncertainty Amplification (by recording site): A term also used for describing the increase xvii in amplitude of seismic waves due to the recording site?s condition. Most seismographs are installed at a site on or near a surface with irregular topography and lithologic structures of heterogeneous material created by erosion, deposition, and other geological processes. These complex near- surface structures tend to amplify the amplitude of incident seismic waves. See also site effect and site response. Annual probability of exceedance: The probability that a given level of seismic hazard (typically some measure of ground motions, e.g., seismic magni- tude, peak ground acceleration or intensity) or seismic risk (typically economic loss or casualties) can be equalled or surpassed within an exposure time of 1 year (see also Hazard curve, Total hazard curve). Area source: See source zone Areal source: See source zone Aseismic: In seismology, an adjective for an area where few or no earthquakes have been observed. In earthquake engineering, an adjective for structures that are designed and built to withstand earthquakes Attenuation relationship: A mathematical expression that relates a ground- motion parameter, such as the peak ground acceleration, to the source and propagation path parameters of an earthquake such as the magnitude, source-to-site distance, or fault type. Its coefficients are usually derived from statistical analysis of earthquake records. Attenuation equation: See also attenuation relationship Azimuth: In seismology, the direction from a seismic source to the seismic Station recording this event. It is usually measured in degrees clockwise xviii from the north. Back-azimuth: The direction from the seismic station toward a seismic source. It is usually measured in degrees clockwise from the north. Base-shear: The horizontal shearing force at the base of the structure. The maximum base-shear, typically in the form of a fraction of the weight of the structure, is an important parameter in earthquake response studies and in earthquake-resistant design (see also seismic hazard zoning coefficient) Converted wave: The conversion of P to S waves and S to P waves that occurs at a discontinuity for non-normal incidence. These converted waves sometimes show distinct arrivals on the seismogram between the P and S arrivals and may be used to determine the location of the discontinuity. Deterministic scenario earthquake: A representation in terms of useful descriptive parameters of an earthquake of specified size postulated to occur at a specified location (typically active fault), and of its effects. Deterministic spectra: It shows the spectral acceleration as a function of spectral period for a specified fractile of the attenuation dispersion. See also deterministic earthquake scenario. Earthquake hazard: See seismic hazard Earthquake magnitude recurrence model: The recurrence rate of earthquake magnitude is assumed to follow the cumulative Gutenberg-Richter frequency-magnitude relationship. Earthquake source parameters: The parameters that are specified for an earth- quake source, depending on the applied model. They are origin time, xix hypocenter location, magnitude, focal mechanism, and moment tensor for a point source model. They include fault geometry, rupture velocity, stress drop, and slip distribution for a finite-fault model Earthquake swarm: A series of earthquakes occurring in a limited area and time period in which there is not a clearly identified main shock of a much larger magnitude than the rest. See also swarm. Epistemic uncertainty: The uncertainty in seismic hazard analysis due to imperfect knowledge in model parameterizations and other limitations of the methods employed. Epistemic uncertainty can be reduced by improvements in the modelling and analysis. See also aleatory uncertainty. Equivalent earthquake magnitude: Magnitude of hypothetical event that would yield energy equivalent to the total energy released in the area Exceedance probability: The probability that a specified value of a strong- motion parameter is exceeded by a future occurrence within a specified period of time. Foreshock: A smaller earthquake preceding the largest earthquake in an earthquake sequence. Hazard curve: It gives the probability that any of the scenario earthquakes will produce a ground motion exceeding the ground motion test value ( See also Total hazard curve) Hawaiian eruption: An eruption of highly fluid (low-viscosity) magma that commonly produces lava fountains and forms thin and widespread lava flows, and generally very minor pyroclastic deposits around the vents. xx Hybrid earthquake: In volcano seismology, an earthquake with characteristics of both high-frequency and long-period earthquakes, often applied to earthquakes that have an impulsive high-frequency onset with onset and S-wave followed by an extended low-frequency coda. INSAR: Interferometric synthetic aperture radar IPGP: Institut de Physique du globe de Paris, France ISC: International Seismological Centre Local site effect: A qualitative or quantitative description of the topography, geology, and soil profile at a site that affect ground motions during an earthquake. Magnitude: Used to define the size of an earthquake. There are many different scales that can be used to define magnitude e.g. moment magnitude (Mw), local magnitude (Ml), short-period body wave magnitude (Mb), surface wave magnitude (Ms). Maximum credible earthquake: The maximum earthquake, compatible with the known tectonic framework, that appears capable of occurring in a given area. Maximum possible earthquake magnitude: Maximum earthquake magnitude that appears capable of occurring in a given area. See also Maximum credible earthquake. Mean Hazard: See seismic hazard analysis and Total hazard NEIC: National Earthquake Information Center in Golden, Colorado, operated by USGS. Paleoseismology: science to study the nature, timing and location of pre- xxi instrumental earthquakes. Peak ground acceleration (PGA): The maximum acceleration amplitude measured (or expected) in a strong-motion accelerogram of an earthquake. Poisson distribution: A probability distribution that characterizes discrete events occurring independently of one another in time. Poisson model: see Poisson distribution Poisson process: A point process in which events are statistically independently and uniformly distributed with respect to the independent variable (e.g. time). See also Poisson distribution. Probability density: A function of a continuous statistical variable whose integral over a given interval gives the probability that the variable will fall within that interval . Probability of Exceedance: The probability that, in a given area or site, an earth- quake ground motion will be greater than a given value during some time period. Receiver function: The spectral ratio of the horizontal component of S waves to the vertical component of P waves recorded at a single station from a teleseismic event. Assuming a horizontally layered structure beneath the station, it gives an estimate of the seismic velocity structure, particularly the nature of discontinuities of the crust and the uppermost mantle. Return period: The average time between exceedance of a specified level of ground motion at a specific location, equal to the inverse of the annual probability of exceedance. xxii Rift (or Rift Valley): An extended crustal feature marked by a fault-caused trough that is created in a zone of divergent crustal deformation. Risk: The probability of harmful consequences, or expected losses (deaths, injuries, damage to property, threats to livelihoods, disruption of economic activity, or damage to the environment) resulting from interactions between natural or human-induced hazards and vulnerable conditions. Conventionally, risk is expressed by the following notation: Risk= Hazard x Vulnerability (Source: www.unisdr.org/eng/library/lib-terminology-engome.html) SAC: Seismic Analysis Code developed by Lee Minner and Peter Goldstein in 1999. Scenario earthquake: See deterministic scenario earthquake SEISAN: A software package for earthquake analysis edited by Jens Havskov and Lars Ottem?ller. Seismic hazard: Any physical phenomena associated with an earthquake (e.g. ground motion, ground failure, liquefaction, and tsunami ) and its effects on land use, man-made structures, and socioeconomic systems that have the potential to produce a loss. It is also used without regard to a loss to indicate the probable level of ground shaking occurring at a given point within a certain period of time. Seismic hazard assessment: The calculation of the seismic hazard, expressed in probabilistic terms, as contrasted with deterministic seismic hazard analysis, for a site or group of sites. The result is usually displayed as a seismic hazard curve or seismic hazard map. xxiii Seismic hazard curve: A plot of probabilistic seismic hazard (usually) specified in terms of annual probability of exceedance) or return period versus a specified ground motion parameter for a given site. Seismic hazard map: A map showing contours of a specified ground-motion parameter or response-spectrum ordinate for a given probabilistic seismic hazard or return period. Seismic source: A localized area or volume generating coherent, usually transient seismic waveforms, such as an earthquake. Seismic zonation: The geographic delineation of areas having different potentials for hazardous effects from future earthquakes. Seismic zonation can be done at any scale: national, regional, local or site. Seismic source zone: An area of seismicity probably sharing a common cause. Seismic zone coefficient: A relative intensity ratio of anticipated earthquake motion taking into account earthquake motion for the entire country or region. Seismic zoning map: A map used to portray seismic hazard or seismic design variables, for example, maps used in building codes to identify areas of uniform seismic design requirements. Seismotectonic source zone: See seismic source zone Seismogenetic source zone: See seismic source zone Site effect: The effect of local geology and topographic conditions at a recording site on ground motion. See also local site effect. Site response: The modification of earthquake ground motion in the time or frequency domain caused by local site condition. xxiv Source zone: An area in which an earthquake is expected to originate; more specifically, an area considered to have a uniform rate of seismicity or a single probability distribution for purposes of a seismic hazard or seismic risk analysis. Swarm: See also earthquake swarm Total hazard: The calculation of the mean annual frequency of exceedance as a function of spectral acceleration for a particular spectral period into a specific site in the study area due to multiple source zones using the total probability theorem. Total hazard curve: It shows the curve of the mean annual frequency of exceedance as a function of spectral acceleration for a particular spectral period at a specific site. Seismic hazard curves are first obtained for each seismic source zone and then combined to obtain the total hazard curve. See also Seismic Hazard curve. Uniform hazard (response) spectrum: In probabilistic hazard analysis, a response spectrum with ordinates having equal probability of being exceeded. Vulnerability: the degree of loss to a given element or set of elements at risk resulting from the occurrence of a specific seismic hazard. It is expressed on a scale from 0 (no damage) to 1 (total loss). Uniform hazard spectra: It shows the spectral acceleration as a function of spectral period for a specified return period. See also Uniform hazard (response) spectrum. xxv LIST OF PERSONAL RECENT ORIGINAL PAPERS This thesis is based on the following original published papers: 1. Mavonga, T., 2007a. Some characteristics of aftershock sequences of major earthquakes from 1994 to 2002 in the Kivu Province, Western Rift Valley of Africa. Tectonophysics, Vol. 439, pp 1-12. 2. Mavonga, T., 2007b. An estimate of the attenuation relationship for the strong ground motion in the Kivu Province, Western Rift Valley of Africa, Phys. Earth. Planet. Interiors, Vol. 162, pp 13- 21. 3. Mavonga, T., Kavotha, K.S., Lukaya, N., Etoy, O., Jacquea, D., 2006. Seismic activity prior to the May 8, 2004 eruption of volcano Nyamuragira, Western Rift Valley of Africa. Journal of Volcanology and Geothermal Research, 158, 355-360. 4. Mavonga, T.G., Durrheim, R.J., 2009. Probabilistic seismic hazard assessment for Democratic Republic of Congo and surrounding areas, Western Rift Valley of Africa, Special AfricaArray issue of South African Journal of Geology, Volume 112, Number 3-4, pp. 329-342 5. Mavonga T., 2009. Crustal structure beneath two seismic broadband stations revealed from teleseismic P wave receiver function analysis in theVirunga Volcanic area, Western Rift Valley of Africa, Special AVCOR issue of Journal of African Earth Science, doi: 10.1016/j.jafrearsci.2009.11.003 6. Mavonga T., Kavotha, K.S., Lukaya, N., Etoy,O., Wafula, M., Rusangiza, K.B., Jacques, D., 2010. Some aspects of seismicity prior to the 27 November 2006 eruption of Nyamuragira volcano and its implication for volcano monitoring and risk mitigation in the Virunga area, Western Rift Valley of Africa, Special AVCOR issue of Journal of African Earth Science, doi: 10.1016/j. jafrearsci. 2010.02.002 1 CHAPTER 1 INTRODUCTION 1.1 General Introduction The Democratic Republic of Congo (DRC) is located between longitude 12?E and 30?E, and latitude 4?N and 12?S. The geology of the DRC is complex. Its most prominent morphological feature is the Western Rift zone that runs along the DRC?s eastern border and through neighbouring countries (e.g. Uganda, Rwanda, Burundi, Tanzania) between 28?E to 32?E and 4?N to 12?S. The Western Rift Valley forms the western branch of the East African Rift System. It extends over 1600 km in an arc-like pattern along Lakes Albert, Edouard, Kivu and Tanganyika, until it joins the eastern branch, after which it extends further to the south along Lakes Rukwa and Malawi. The rift is generally bounded by tilted blocks of crystalline Precambrian rocks [Morley, 1988] that are cut by tertiary faults that strike parallel to the rift axis. The Western Rift Valley of Africa (WRA) has experienced several severe earthquakes and volcanic eruptions in recent historical times. The earthquakes with magnitude M ? 6 occurred mostly in DRC and neighbouring countries (e.g. Uganda and Tanzania). They were located in the northern Lake Tanganyika region (e.g. 22 September 1960, Ml6.5), the Ruwenzori area (e.g. 20 March 1966, Ml7.0; 5 February 1994, Ms6.5), the region southwest of Lake Tanganyika (Kabalo earthquake, 11 September 1992, Ms6.7), and the Central Kivu basin (Kalehe earthquake, 24 October, 2002, Mw6.2). In general, seismic activity in the Western Rift Valley is confined to the following zones: the south, central and 2 northern part of the Lake Tanganyika, the western border of Lake Kivu and Ruzizi Valley, the Lake Edouard trough, and the Mount Ruwenzori area (Figure 1.1). Furthermore, temporary seismic observations made in an area west of the Rift Valley near Lake Kivu show that seismic activity is not only high in the Rift Valley, but also outside the Rift Valley, especially in the Masisi area located about 50 km west of the Western Rift Valley [Zana et al., 1990]. On 29 April 1995, an earthquake of magnitude Mb5.1 occurred in the Masisi area. Earthquakes with magnitudes greater than 5 were also observed in the Congo intracratonic basin, located about 450 km from the rift axis [Wafula and Zana, 1990; Zana et al. 1992 and Atalay, 2002] and adjacent areas (northeast region of Republic of Congo Brazzaville). Seismic activity in the basin of Lake Kivu was moderate (Magnitude ? 4) until 1997, except for the occurrence of an earthquake of magnitude Ml4.5 on 25 April 1965 with its epicentre 10 km south of Bukavu City [Wohlenberg, 1968] and another on 6 January 1976 (Mb5.2) in the Ngweshe area. Since 1997, landslides have started to occur in Bukavu City and vicinity due mainly to seismic activity observed in the southern part of Lake Kivu (e.g. 11 February 1997, Mb4.7; 11 November 1997, Mb4.7; 2 March 2000, Mb5.4; 3 March 2000, Mb4.9; 12 September 2000, Mb4.6; Munyololo et al., [1999]). 3 Figure 1.1 Large earthquakes with magnitude M?6 observed in the Western Rift Valley of Africa since 1910. The coordinates of the 1910 M7.3 were obtained by combining those obtained from instrument with local macro seismic record as published in Ambraseys and Adams [1992] LA= Lake Albert, R= Ruwenzori Mountain, E= Lake Edouard, K=Lake Kivu, TA= Lake Tanganyika, LM= Lake Moero and RV= Ruzizi Valley On 24 October 2002, an earthquake of magnitude Mw6.2 occurred in the central part of Lake Kivu. Aftershocks of this earthquake extended from the Kalehe area to Idjwi Island. Two people were killed and there was significant damage to 4 masonry buildings [Mavonga, 2007a]. Recently, on 3 February 2008, an earthquake of magnitude Mw6.0 occurred about 20 km north of Bukavu City. Nine people were killed, more than one thousand houses destroyed and more than 400 people were injured (Figure 1.2). Figure 1.2 Damage to a house resulting from earthquake on 3/2/2008 at Bukavu City (D.R Congo). Source: I. Badriyo, Project UNOPS ?Unite de gestion de risques volcaniques?, Goma, North-Kivu, Democratic Republic of Congo Furthermore, some earthquakes have also occurred in the neighbouring countries (Rwanda and Burundi). On 20 March 2003 and 24 February 2004, earthquakes of magnitude Mb5.2 and Mb4.8 occurred to the southeast of Lake Kivu in Kibuye territory in Rwanda and about 35 km from Bujumbura in Burundi, respectively. Two people were killed by the latter earthquake. 5 The overwhelming majority of the shocks recorded since 2002 have occurred on normal faults of the western Kivu border fault system. It is noteworthy that in the northern edge of Lake Kivu, inhabitants of Goma City and surrounding areas (e.g. Gisenyi in Rwanda) are also threaten by volcanic hazard due to eruptions of the Virunga volcanoes (Nyiragongo and Nyamuragira). 1.2 Motivation for this study The seismic hazard in the Democratic Republic of Congo and areas adjacent to the Western Rift Valley of Africa has been previously assessed by Midzi et al. [1999], Twesigomwe [1997] and Zana et al. [1992a]. Due to the large area of approximately 50? x 25? that was investigated, Midzi et al. [1999] considered only regional structures in their seismic source zonation and did not take the local details of observed seismicity and tectonic features in the Western Rift into account. The Twesigomwe [1997] study was mostly based on the assessment of seismic hazard in Uganda. The estimation of earthquake hazard in Zaire (former name of DRC) by Zana et al. [1992a] was based only on combining the pattern of spatial distribution of epicentres, the equivalent earthquake magnitude distribution and trends of a and b in the Gutenberg Richter formula, but they did not take the attenuation of ground acceleration into account. The seismic-volcanic activity in the Virunga volcanic area has been monitored since 1983 by the Goma Volcano Observatory (GVO). The observation network 6 of GVO is composed of stations Katale (KTL), Luboga (LBG), Kunene (KNN), Rusayo (RSY), Kibumba (KBB), Goma and Bulengo (BLG). The determination of the earthquake hypocenter is a powerful means to locate magmatic activity in the Virunga area. However, locations may be biased due to the lack of a P- and S- wave velocity structure specific for the Virunga area [Chiarabba et al., 2000; Zhang and Thurber, 2003; Stephen, 2005]. Data on deeper structure in the Western Rift are extremely sparse [e.g. Prodehl and Mechie, 1991]. There are no deep seismic reflection or refraction data available for the Western Rift near the Virunga volcanic area. Thus, it is imperative to image P- and S-wave velocity structure beneath Virunga and adjacent areas. The social conditions in the Great Lakes region are continually changing and planning based on seismic and/or volcanic hazard consideration is poor to non- existent. Population growth has led to relatively uncontrolled use of land for building and the development of unsuitable sites that are vulnerable to earthquake and/or volcanic hazard. For these reasons, a revised assessment of hazard is urgently needed everywhere within these areas in DRC and surrounding countries. This will provide us a means of minimizing the danger to life as population density continue to rise. This work will contribute to: 1. The formulation of a national seismic design code for buildings that will assist architects and engineers to take the seismic hazard into consideration, and 7 2. An understanding of how volcanoes work and thereby reduce the risk posed by volcanic eruptions in the Virunga area by deriving a new crustal seismic velocity model for the Virunga area in order to map the migration of hypocenters of earthquakes accurately and identified trends that could be precursors of volcanic eruptions 1.3. Aim and structure of the study The major output of this study is a seismic hazard map for the DRC and adjacent areas based on probabilistic hazard computation using the instrumental earthquake data for a period of about 90 years. The basis for all seismic hazard assessment is the analysis of seismicity or the occurrence of earthquakes in space and time. Accurate knowledge of seismicity is an important tool for understanding active tectonics. One of the basic elements in assessing seismic hazard is the identification of seismic source zones that could affect the particular location at which the hazard is being evaluated. These seismic source zones are also called seismogenetic or seismotectonic source zones. Defining and understanding seismotectonic source zones is often the major part of a seismic hazard analysis and requires knowledge of the regional and local geology, seismicity and tectonics. In this thesis, due to the lack of well documented knowledge of active faults or linear sources in the study area, we used seismotectonic source zones which appear as area sources. 8 Chapter 2 is a review of previous studies on the tectonics and seismicity of the study area. We divide the study area into three seismic source zones: 1. The Western Rift Valley of Africa 2. The Upemba and Moero Rift, and 3. The Congo basin Based on local geology and seismicity, the Western Rift Valley was subdivided into 4 sub-zones: 1. Southern Sudan, Ruwenzori and Lake Edouard trough 2. Virunga volcanic complex, Rutsuru and Masisi area 3. Basin of Lake Kivu, Ngweshe and Ruzizi plain, and 4. Lake Tanganyika Rift In chapter 3, the seismic characteristics of the study area are modelled as a Poisson process following the standard engineering seismology assumptions. The parameters used to characterise each seismic source zone are: 1. Average rate of occurrence or mean seismic activity rate ? (which is a parameter of the Poisson distribution), 2. Level of completeness of the earthquake catalogue (Mmin), 3. Maximum possible earthquake magnitude (Mmax), derived using the method described by Kijko [2004], 4. Gutenberg-Richter [1954] ?b-value? (which indicates the relative number of large and small earthquakes, ?=b ln10), 5. Focal depth, and 6. Regional attenuation relationship for the strong ground motion. 9 These parameters were determined and used as input parameter for hazard computation. The EZ-Frisk package (Version 7.24) [Risk Engineering Inc, 2007] was used to compute the seismic hazard at sites of interest for various combinations of seismic source zone, attenuation equations and two focal depths. Results were displayed in chart of the total hazard curve including the contribution of hazard by seismic source. The PGA was calculated for a specified return period or annual chance of exceedance, on a grid covering the study area specified by their longitude and latitude at a grid interval of 0.5 degree. These values were used to create a seismic hazard map. Maps show the PGA (in unit g) for 2%, 5% and 10% probability of exceedance in 50 years which corresponds to return periods of 2475, 975 and 475 years, respectively. The highest seismic hazard was obtained in the Lake Tanganyika Rift where the PGA values of 0.32g, 0.22g, and 0.16g are expected to be exceeded with probability 2%, 5% and 10% in 50 years, respectively. We then discuss the source of uncertainties in the input parameters used for hazard computation and assess the sensitivity of the outputs to these uncertainties. We also evaluate the effect of delimitation of seismic source zones and the computation method (probabilistic or deterministic hazard investigation) on the final result of seismic hazard analysis. Finally, we compare our results with recent work on earthquake hazard assessment in the Western Rift Valley. 10 In chapter 4, in order to understand how the Virunga volcanoes work and reduce the risk due to volcanic eruptions, we attempt to improve the seismic velocity model beneath the Virunga area. Teleseismic P-waveform receiver function analysis was used to study the crustal structure beneath the area sampled by two broadband stations in the Virunga area. The new velocity model obtained was used to recalculate the hypocenters for the period 18 August 2002 to 7 May 2004 prior to the 8 May 2004 Nyamuragira eruption and 1 July 2004 to 27 November 2006 prior to the 27 November 2006 Nyamuragira eruption. In chapter 5, we synthesize the research findings: a) The seismic hazard assessment for DRC and surrounding areas was used to identify four hazard zones as follows: Zone A (very high hazard), Zone B (high hazard), Zone C (moderate hazard) and Zone D (low hazard). Zone A includes the Lake Tanganyika Rift. b) We indicate how the improved hypocenter locations in the Virunga area, integrated with other geophysical, geological or geochemistry data, may be used to characterize volcanic process and forecast volcanic eruptions. Some recommendations were also made for the future studies on earthquake hazard assessment in the study area and crustal structure beneath Virunga area. 11 CHAPTER 2 GEOLOGY AND SEISMOTECTONICS OF THE DRC AND ADJACENT AREAS 2.1 Introduction The occurrence of earthquakes in DRC and adjacent areas is mainly controlled by the Western Rift Valley of Africa (WRA). A concentration of epicentres follows the rift structures between southern Sudan and southern Malawi. The WRA is occupied by the Lakes Albert, Edouard, Kivu, Tanganyika, Rukwa and Malawi. Furthermore, the region south of the Lake Tanganyika Rift, which includes the Katanga province of south-eastern DRC and north-western Zambia, is of considerable tectonic interest since geological and geophysical studies have revealed seismically active areas that may be related to the WRA with a common NE-SW trend. Although features related to active faults are poorly exposed at the surface, the region is marked by a large negative Bouger anomaly [Mondeguer et al., 1989; Fairhead and Girdler, 1972; Fairhead and Henderson, 1977; Sebagenzi et al., 1993; Shudofsky, 1985]. The following young rifts have been observed in that region: the NE-SW trending Upemba, Luano, Lukusashi, Luangwa Karoo and Moero Cenozoic Rifts. Diffuse seismicity is also observed in the Congo basin [Atalay, 2002]. Studies of local seismicity from micro-earthquake networks and seismotectonics in the Western Rift Valley of Africa, Congo basin, south-eastern DRC and northwestern Zambia have been carried out by De Bremaecker [1955,1959], Wohlenberg [1968], Bram [1972], Maasha [1975], Zana [1977], Zana and 12 Hamaguchi [1978], Zana et al. [1989], Zana et al. [2004], Maasha and Molnar [1972], Tanaka et al. [1980], Zana and Tanaka [1981], Shudofsky [1985], Wafula and Zana [1990], Foster and Jackson [1998], Atalay [2002], Sebagenzi and Kaputo [2002] and Mavonga [2007a]. These studies, either of the focal mechanism of individual events or composite focal mechanism solutions of micro-earthquakes, have shown that the stress field is predominantly of normal faulting type throughout the region, with the exception of the Congo basin. 2.2 Main seismic source zones for the study area The main seismic source zones currently known for the study area are shown in Figure 2.1 to 2.7 and summarized as follows: 2.2.1 The Western Rift Valley of Africa (WRA) A hypsographic Digital Elevation Model (DEM) of the East African Rift System (EARS) is shown in Figure 2.1. Based on local seismicity and geological structure, four seismic sub-zones have been identified in the WRA. (i) Southern Sudan, Ruwenzori area, and Lake Edouard trough (Figure 2.2), (ii) Virunga volcanic complex, Rutsuru and Masisi area (Figure 2.3), (iii) Lake Kivu basin, Ngweshe area and Ruzizi plain (Figure 2.4), and (iv) Lake Tanganyika Rift (Figure 2.5). The main features of each sub-zone are described below. 13 Figure.2.1 Hypsographic DEM of the East African Rift system. Black lines are faults. Source: Chorowicz, J. [2005]. 14 2.2.1.1 Southern Sudan, Ruwenzori area, and Lake Edouard trough The South Sudan is dominated by relatively strong earthquakes (e.g. sequence of large events in 1990-91) but with poor tectonic control [Girdler and McConnell, 1994]. The Ruwenzori mountain is an uplifted horst of Pre-Cambrian rock [Cahen, 1954] bounded to the north by the NNE-SSW trending Lake Albert basin and to the south by Lake Edouard [Ebinger, 1989b] (Figure 2.2). Figure 2.2 Generalised geology map of the Kivu Province. AFZ=Aswa fault; LK=Lake Kyoga; LA=Lake Albert; R=Ruwenzori Mountain; KB=Katonga Break; RFB=Ruwenzori Fold Belt; E=Lake Edouard; MA=Masisi area; V=Virunga Volcanic Field; K= Lake Kivu; SK=South Kivu Volcanic Field; T=Lake Tanganyika.Source: Modified from Twesigomwe [1997] and Mavonga [2007b]. 15 The Ruwenzori area experienced large earthquakes on 20 March 1966 (Mw=6.8) and 5 February 1994 (Mw=6.2) that killed 160 and 8 people, respectively [National Earthquake Disaster Committee, 1994; Mavonga, 2007a]. Lake Edouard is also an active area. Many felt earthquakes have been reported by local inhabitants at Butembo City (DRC) located at the northern end of Lake Edouard. The last important earthquake (Mb=5.2) in the Lake Edouard area occurred on 8 May 2003. 2.2.1.2 Virunga volcanic complex, Rutsuru and Masisi area The Virunga volcanic complex is the largest of the Cenozoic volcanic complexes in the Kivu Province and the only one that is presently active. The Virunga volcanics are divided into three subgroups: the Eastern (Muhavura, Gahinga, Sabinyo), Central (Visoke, Karisimbi, Mikeno), and Western (Nyiragongo and Nyamuragira). Except for a short lived eruption of Visoke in 1956, only the Western group has had eruptive episodes in recent historical times [Mavonga et al., 2006, Kavotha et al., 2003, Komorowsky et al., 2003]. The Virunga volcanoes are distributed along an east-west trend at right angles to the rift axis [Kampunzu et al., 1986]. The tectonic seismicity in the volcanic area is very low (M? 4) (Figure 2.3). The Rutsuru basin is located at the northern part of the Virunga area (Figure 2.3). It is bordered by fault segment. Felt earthquakes are often reported by the local inhabitants in the Rutsuru area. 16 Masisi is located at the northwest of Lake Kivu (Figures 2.3 and 2.4). A study of earthquake focal mechanisms by Tanaka et al. [1980] showed that the direction of the fault traces in that area is SE-NW, and the average focal mechanism is normal faulting with the tension axis perpendicular to the strike of the fault traces. The last strong earthquake occurred in the Masisi area on 29 April 1995 (Mb=5.1) [Mavonga, 2007a]. 2.2.1.3 Lake Kivu Basin, Ngweshe area and Ruzizi plain The Kivu basin consists of two subsiding half grabens separated by the 700 m high Idjwi horst structure (Figure 2.3). The basin is bordered to the west by the prominent (more than1500 m high) escarpment of the West Kivu border fault segment, and to the east by the prominent escarpment of the East Kivu border fault segment [Ebinger, 1989a]. A network of oblique-slip transfer faults trending transversely to the main rift axis accommodates differences in elevation along the rift axis, and has contributed to the very complex tectonic relationship. Moreover, the Lake Kivu basin marks the transition of the Western Rift Valley from the predominant NW-SE orientation of Ruzizi-Tanganyika Rift basin to the more SW-NE trend of the Kivu, Virunga and Rutsuru-Lake Edouard Rift basins [Ebinger, 1989b]. The central part of Lake Kivu experienced a large earthquake (Mw=6.2) on 24 October 2002 in the Kalehe area, which was felt strongly at Goma, Bukavu and Kigali. This earthquake was the largest observed in the basin of Lake Kivu since 1900 [Mavonga, 2007a]. On 03 February 2008 an Mw=6.0 earthquake occurred 17 20 km north of Bukavu City. This earthquake claimed 9 lives, more than 400 people were seriously injured and about 1500 houses collapsed. Figure 2.3 Map showing the eight volcanoes of the Virunga Volcanic Field and seismicity of Lake Kivu and surrounding area for the period 1976-2002. (NY= Nyiragongo; NM= Nyamuragira; KA= Karisimbi; MI= Mikeno; VI= Visoke; SA= Sabinyo; GA= Gahinga; MU= Muhavura). The filled triangles and squares indicate seismographic stations and earthquake epicentres, respectively. The filled small circles indicate cities [Modified from Komorowsky et al., 2003 and Mavonga, 2007a] 18 The Ngweshe area is located to the southwest of Lake Kivu (Figure 2.4). Its physiography is characterized by dislocated highlands of Precambrian age without indication of recent volcanism. Several felt earthquakes have been reported in the Ngweshe area. The trend of tension axis in this area is SE-NW perpendicular to the fault trace [Tanaka et al., 1980]. The Ruzizi plain is located between Lake Kivu and north of Lake Tanganyika. Three strong earthquakes occurred in the Ruzizi area (M(Lwiro)=5.5 in 1954, 5.0 in 1956, 6-6.5 in 1960). Figure 2.4 Epicenter map for the period 1955 to 1962 showing the central part of the Western Rift Valley of Africa (Masisi, Lake Kivu, Ngweshe and Ruzizi plain). The tectonic lines are from Cahen [1952] and Furon and Dumain [1959]. The filled circles, rectangles and triangles are earthquakes with M?5, 4?M<5 and 3?M<4, respectively. Source: J. Wohlenberg [1968]. 19 2.2.1.4 Lake Tanganyika Rift Many larger magnitude earthquakes have occurred in the Tanganyika Rift (Figure 2.5). The best known are: The event of 13 December 1910 at the southern end of Lake Tanganyika (M=7.3), The 22 September 1960 Uvira earthquake (Ms=6.5) [Zana and Hamaguchi, 1978], with epicentre at the northern end of Lake Tanganyika, and The 5 December 2005 earthquake (Ms=6.8), which took place in the central part of Lake Tanganyika. It was felt even in Nairobi (Kenya). The accurate history of the seismic activity of this rift is not known. However the topography of this area is marked by grabens and horsts that are most likely associated with earthquakes [Mortelmans, 1951; Zana et al., 2004]. 2.2.2 Southeastern DRC and northwestern Zambia The most prominent seismotectonic features in this region are the Upemba and Moero or Mweru Rifts (Figure 2.5). The Upemba Rift is characterized by a NE- SW striking fault extending along its eastern side [Studt et al., 1908]. The Upemba Rift may extend northward to the Kabalo area, which experienced an earthquake with magnitude Mw6.5 on 11 September 1992. The earthquake parameters of the main shock provided by USGS are: origin time 03:57:26.2(UT), latitude 06.091?S, longitude 26.680?E, depth 10 km (fixed). In the city of Kabalo, the poorly constructed brick buildings experienced the most severe damage. The 20 main shock claimed 11 lives and 109 people were seriously injured (4 died within a month of the main shock). More than 2000 families were left homeless. Detailed investigation has revealed that the main geological features in the Kabalo area trend in the NNE-SSW direction, similar to those found in the Upemba Rift [Zana et al., 2004]. Sebagenzi and Kaputo [2002] reviewed the available gravity, heat flow and seismological data in the region, and found geophysical evidence of continental break-up in the southeastern DRC and northwestern Zambia. Continental rifts preferentially follow pre-existing basement weaknesses such as ancient orogenic belts [e.g. Dunbar and Sawyer, 1989; Vauchez et al., 1997]. In the southeast of the DRC and northwestern Zambia these rifts are associated with belts which developed during the Paleo- and Mesoproterozoic times. The following rifts are known: the NE-SW trending Upemba, Luano, Lukusashi and Luangwa and Moero rifts. The stratigraphy of the Luangwa rift basin is relatively well known, but that of the Moero rift is practically unknown. In the Luano rift, Late Carboniferous- Permian strata appear to be present, but the presence of Early Jurassic strata has still to be demonstrated [Delvaux, 2001]. Little is known about the Upemba rift, which is located entirely in the DRC. At the present state of knowledge, only some Upper Carboniferous sedimentary rocks succeeding a basal tillite and Cenozoic terranes are known to have filled the Upemba rift [Mortelmans, 1953; Villeneuve, 1983]. The opening of the NE-SW trending Upemba and Luano- Luangwa-Lukusashi Karoo basins was related to crustal extension during the breakup of the Gondwana [Daly et al., 1989], while the Moero was opened during 21 the Cenozoic rifting episode [Mondeguer et al., 1989] (See Figure 2.5 and Figure 2.6). According to seismological data, both Karoo and Cenozoic Rifts, which are still opening, may correspond to incipient arms of the East Africa Rift System (EARS) [e.g. Camelbeeck and Iranga, 1996; Fairhead and Girdler, 1972; Fairhead and Henderson, 1977; Lombe and Mubu, 1992]. This crustal extension has a tensional stress axis perpendicular to the NE-SW trending axis of the regional gravity low parallel to the lithosphere thinning axis [e.g. Fairhead and Henderson, 1977; Sebagenzi et al., 1997; Sebagenzi et al., 1993] (Figure 2.6). The b-value of the Upemba and Moero Rifts seismic source zone determined from the seismic history of the area in this study (chapter 3) is identical, within error limits uncertainty, with that of the entire Western Rift Valley seismic source zone from south of Sudan-Ruwenzori to the Lake Tanganyika Rift. Also, Mavonga [2007a] showed a correlation between fault area and magnitude of earthquakes which occurred in the Western Rift Valley and that of the Upemba area. Therefore, the evidence suggests that Upemba and Moero Rifts are part of the Western Rift Valley, but additional data are needed to confirm this as the relationship between the surface structures and deeper features is not known [Sebagenzi and Kaputo, 2002]. 22 Figure 2.5 A map of Eastern African Rift System (EARS) showing the continental breakup in the southeast of DR Congo and northern Zambia (Modified from Kinabo et al., 2007) The small red stars indicate large earthquake with magnitude M?6 observed in the Western Rift Valley of Africa since 1910. 23 Figure 2.6 Regional gravity anomaly and heat flow determination values at the sites reported from Zambia (solid triangles) by Chapman and Pollack[1977], Kenya and Tanzania (solid squares) by Nyblade et al. [1990] and from D.R. Congo (encircled stars symbols) by Sebagenzi et al. [1993]. Barbed lines represent active normal faults; labeled thin lines and unlabeled dashed lines represent gravity anomalies with a contour interval of 5 mgal. [Modified from Sebagenzi and Kaputo, 2002] 24 2.2.3 The Congo basin Due to the scarcity of outcrops and a very dense equatorial forest coupled with a thin cover of recent unconsolidated sediments over most of the area, the geology of this basin is studied in a few localized areas (e.g. Daly et al. 1992). No surface fault displacements have been documented, even though some large and damaging shocks have occurred. Studies of four earthquakes with magnitude ranging from Mb5.4 to Mb5.6 that occurred in the Congo basin during the period 1976 to 1998 (Fairhead and Stuart [1985], Diewonski et al. [1996], Atalay [2002]), together with previous investigations, demonstrate that the Congo basin is in a state of predominantly horizontal compression. The fault mechanisms of these earthquakes show approximately E-W oriented P-axes, which could be interpreted as E-W contraction of the African Plate due to ridge push forces from the Mid-Atlantic Ridge and the East African Rift System (EARS) (Figure 2.7) 25 (a) (b) Figure 2.7 (a) Map showing the regional setting of the Congo Basin (b) Focal mechanism plots of all studied earthquakes for source dynamics in the Congo Basin. The black circles are the reported seismicity of the area by USGS in its entire history of reporting period whilst the white overlapping dots are the earthquakes studies by Atalay (2002). The fault plane solution for the 980305 earthquake is done from P-wave first-motion data and the mechanism for the 980426 event is further constrained by P- and SH-waveform fitting. [After Atalay, 2002] 26 CHAPTER 3 EARTHQUAKE HAZARD ASSESSMENT FOR DRC AND SURROUNDING AREAS 3.1 Data All seismic data used in this study are instrumental data compiled from various sources covering the region 14o S to 6o N and 10o E to 32o E for the period 1910 to 2008. The area of this study includes almost all Democratic Republic of Congo excluding the western region, as well as parts of Uganda, Tanzania, Rwanda, Burundi, Sudan and Zambia. Reports of earthquakes were searched for carefully as far back in time as possible using all available catalogues and seismological bulletins published by various reporting agencies. The main sources of data for earthquakes after 1953 are the seismological bulletins published by the ?Institut pour la recherche en Afrique Centrale (IRSAC)? for the period 1953 to 1977, and by its successor, the ?Centre de Recherche en Sciences Naturelles (CRSN)? for the period 1977 to 2007. The seismograph network of the IRSAC in Eastern DRC has been in operation since May 1953 when the station at Lwiro was set up. This network operated initially with three stations: Lwiro (LWI), Butare (BTR) and Uvira (UVI) [Sutton and Berg, 1958]. It was extended later by setting up the Rumangabo (RMG) and Butembo (BTC) stations in the North Kivu Province. In cooperation with the IRSAC, the ?Union miniere du Haut Katanga (UMHK)?, now ?Generale de carrieres et des mines (GECAMINES)? set up a seismological network for Katanga from 1960 to 1970. This network included the stations of Delcommune 27 (DCC), Lubudi (LBD) and Mulungwishi (MLN) (Bram, 1972). However, most of these seismic stations were closed during the period 1964 to 1970 due to the political instability in DRC. Current seismograph stations in DRC are operated by the ?Centre de recherche en Sciences Naturelles (CRSN)? and the University of Lubumbashi (UNILU). At present there are 9 stations, mostly concentrated in the Virunga volcanic area. In July 2007, the ?Musee Royal de l?Afrique Centrale (MRAC)? provided a broadband seismic station at Lubumbashi, in collaboration with UNILU and AfricaArray. Apart the IRSAC/CRSN network, data were obtained from the catalogues of the following reporting agencies: The International Seismological Centre (ISC, previously International Seismological Summary, ISS), United States Geological Survey (USGS, which includes NEIS, NEIC, USCGS and CGS), Zimbabwe Meteorological Service Seismological Bulletin, Bulawayo (BUL), Gutenberg-Richter [1949, 1954] and Ambraseys and Adams [1992] 28 3.2 Method 3.2.1 Unification of magnitudes A catalogue was compiled in Seisan format, listing information on the source of the data, the date, origin time, coordinates and magnitude of the earthquake. For source parameters other than magnitude, the most reliable and likely solution was selected if more than one solution was reported. The choice was made in order: ISC, USGS, LWI and BUL. However, there are several different magnitudes scales in use today. Even when magnitudes are determined from the same scale, there often are significant differences in the magnitude determined by different seismological centres. The reasons for the difference lie in the definition of magnitude (B?th, 1981). It is therefore, necessary to homogenize the magnitude by choosing a reliable magnitude scale before carrying out any seismic study. We decided to use the moment magnitude Mw for our unified catalogue, as the moment magnitude Mw is a direct indicator of the seismic moment of an event [Boore and Joyner, 1984; Joyner, 1984]. For simplicity, the USGS and BUL magnitudes were converted to the ISC magnitude using the least squares regression relations of Hlatywayo (1992, 1995): Mb (ISC) = 0.90?0.14 Mb (USGS) +0.38?0.06 Mb (ISC) = 0.59?0.08 Mb (BUL) +1.97?0.26 (3.1) Both the body wave and surface wave magnitudes have been reported by ISC. From this catalogue, a relationship was obtained between Ms(ISC) and Mb(ISC) as shown in Figure 3.1 : Ms(ISC) = 1.358 Mb(ISC) -2.311 r2 = 0.56 (3.2) 29 Figure 3.1 Regression of Ms versus Mb ISC?s magnitudes for events observed in the Western Rift Valley of Africa The magnitudes published by Lwiro (LWI) station were determined by a procedure similar to that proposed by Richter [1935] using the records of standard torsion seismometer by Anderson and Wood [1925]. The Nordquist [1945] nomogram, adapted for the short period Benioff seismometer (To=1 sec, Tg=0.25 sec, Mag=100K, [De Bremaecher, 1955]), was used to determine LWI magnitudes. For the period from 1953 to 1964, the Lwiro magnitude M(LWI) reported by the IRSAC catalogue was the only data available for the unified catalogue. For magnitudes less than 4, M(LWI) is more in line with Richter local magnitude. However, for magnitudes greater than 4, M(LWI) shows a big departure from Richter local magnitude. In this range, M(LWI) gives magnitudes which are comparable to teleseismic body wave magnitude Mb [B?th, 1975]. So, to obtain a reliable magnitude scale, M(LWI) was converted to an equivalent USGS Mb 30 magnitude. Using 86 selected earthquakes of magnitude ?4 for the period from 1965 to 1977 recorded by the short-period Benioff seismograph at Lwiro seismological station and those provided by USGS (Table A.1), a regression relationship was determined between Lwiro magnitude M(LWI) and USGS?s Mb magnitude using the least square method as: Mb (USG) = 3.315+0.282 M(LWI) r2 =0.521 (3.3) The corresponding regression curve is shown in Figure 3.2. Figure 3.2 Regression of M(LWI) versus Mb (USGS) It is noteworthy that, statistically, the linear regression lines obtained by equations (3.2) and (3.3) are poor due to limited data available used for this regression analysis. The linear relationship gets stronger when the coefficient of determination r2 gets closer to 1. 31 The moment magnitude Mw was chosen as a unified magnitude in the compiled catalogue. The following relations were used to convert various magnitudes to Mw. For small events, local magnitude Ml is considered to give a reliable measure of event size [Singh et al., 1990]. Hanks and Kanamori [1979] showed that surface wave magnitude Ms, when converted to magnitude Mw, has a similar relation to Ml. The relations used to convert Ml and Ms to moment magnitude Mw via seismic moment Mo are: Log10 Mo =1.5Ms +(16.1?0.1) Log10 Mo =1.5Ml +16.0 (3.4) where, Mo is seismic moment given in dyne-cm and moment magnitude Mw is given by the relation: Mw = 2/3 log10 Mo -10.7 (3.5) The conversion relation is valid for the range 3? Ml ? 7 and 5? Ms ? 7.5. To convert Mb to Ms a theoretical relationship was adopted, the Ms-Mb relation [Marshall, 1970]: Ms = 2.08 Mb ?5.65 (3.6) To convert Ml to Mb, we used the relationship of Richter [1958]: Mb=1.7 +0.8 Ml ?0.01(Ml) 2 (3.7) 32 3.2.2 Completeness of the catalogue The statistical distribution of a magnitude for a group of earthquakes follows the Gutenberg-Richter [1949] frequency-magnitude relationship: Log10 N(M)= a- bM ; (3.8) where N is the number of earthquakes of magnitude M or greater, and a and b are constants. In this equation, a is a measure of seismicity dependent of the area and time involved and b is a measure of the relative abundance of large and small earthquakes. The b-value may vary between 0.7 to 1.1 [Frohlich and Davis, 1993; McGuire, 1993]. The completeness of a catalogue is defined as that range over which the magnitude-frequency relation roughly follows the Gutenberg-Richter relation. The highest magnitude (at the low magnitude end) that departs from the linear range is called the threshold magnitude. This discrepancy at the low magnitude end is attributed to incompleteness of detection. The most widely used approach to determine the threshold magnitude is to apply the cumulative frequency distribution instead of single frequency distribution. This method also smoothes out error in the magnitude [Hlatywayo, 1992]. The magnitude-frequency distribution of the 2249 earthquakes in this catalogue is shown in Figure 3.3. Assuming the linearity of the magnitude relationship above magnitude 3.7, the plot implies that the threshold of completeness is around magnitude 4.0 . 33 Figure 3.3 Frequency distribution of M?2 earthquakes in the 1910-2007 catalogue. Magnitude increment is (a) ?m=0.5 and (b) ?m=0.25. The filled squares and lozenges are number of events and accumulated number of events, respectively. 34 3.2.3 Probabilistic seismic hazard analysis In this study, the probabilistic approach to seismic hazard analysis as formulated by Cornell [1968] and McGuire [1976, 1993] was employed. Application of the procedure includes several steps: 1. Potential seismic source zones are defined, over which all available information may be averaged. These zones are usually associated with geological or tectonic features (e.g. faults). 2. Seismicity parameters are determined for each seismic source zone. Assessment of the above parameters requires a seismic event catalogue containing origin times, size of seismic events and spatial locations. Site- specific analysis of seismic hazard requires knowledge of the attenuation of the selected ground motion parameter, usually PGA, as a function of earthquake magnitude and distance. 3. The final step requires the integration of individual contributions from each seismic source zone into a site. A Poisson model of earthquake occurrence, which assumes that events are independent, was adopted [Bender and Perkins, 1987]. Therefore, foreshocks, aftershocks, and earthquake swarms were removed from the initial catalogue of 2249 events. Furthermore, Mw=4.0 was selected as the lower magnitude bound (Mmin) because smaller earthquakes are considered unlikely to cause damage, even to houses that are poorly designed and built (as is typical for the DRC and adjacent areas). Thus any remaining events with M<4 were also excluded from the catalogue, leaving a sub-catalogue of 821 events. Scientists wishing to use the data should direct a request to CRSN/Goma Volcanic Observatory. 35 3.2.4 Seismic source zones Although most earthquakes are caused by faulting, it is not possible to identify individual seismically-active faults owing to the lack of good coverage of seismographic stations in the region and/or detailed Quaternary fault maps at a suitable scale showing the active faults. The characteristic elements, such as fault geometry, slip rate of faulting during recent geologic period, and fault segmentation length are not available. However, our current knowledge of the tectonics and seismicity of the study area (chapter 2) is good enough to designate source area models within which earthquake characteristics may be assumed to be uniform. These previous seismotectonic studies suggest that three main seismic source zones are responsible for the damaging earthquakes in DRC and surrounding areas, namely: 1. Upemba-Moero Rifts (an area encompassing the south-eastern part of DRC and north-western Zambia), 2. Congo Basin, and 3. Western Rift Valley of Africa, which was subdivided into 4 sub-zones on basis of the local seismicity and geological structure (i) Southern Sudan, Ruwenzori area and Lake Edouard trough, (ii) Virunga volcanic complex, Rutsuru and Masisi area, (iii) Lake Kivu Basin, Ngweshe and Ruzizi plain, and (iv) Lake Tanganyika Rift. See Figure 3.4. 36 (a) (b) Figure 3.4 (a) Epicentres in the D. R. Congo and surrounding areas for the time period 1910-2007 used in the seismic hazard analysis. The filled large stars, rectangles, lozenges and small stars indicate earthquakes of magnitude range between 4.4 to 7.3, 4.3 to 4.4, 4.2 to 4.3 and 4.0 to 4.2, respectively. (b) The seismic source areas are marked as polygons (1=Upemba-Moero Rift; 2=Congo Basin; 3=Western Rift Valley). The Western Rift Valley seismic source zone is divided into four sub-zones (3a=Ruwenzori-Lake Edouard trough; 3b=Virunga volcanic complex-Rutsuru-Masisi; 3c= Lake Kivu Basin-Ngweshe-Ruzizi; 3d=Lake Tanganyika Rift). 37 All these sub-zones are assumed to be governed by the same tectonic setting, although the southern Sudan may not be actually part of the Rift, but it seems similar to others in the sub-zone (i) in terms of seismicity. The geographic boundaries of these seismic source zones and sub-zones are shown in Figure 3.4, superimposed on seismicity. 3.2.5 Seismic parameters of seismic source zone Following typical assumptions made in engineering seismology, the seismic characteristics of each seismic source zone were modelled as a Poisson process. The most widely used model of magnitude distribution has its source in the classical Gutenberg-Richter relation. In order to ensure a finite seismic energy release [Knopoff and Kagan, 1977], the Gutenberg-Richter relation is often combined with an assumption on the existence of a physical upper limit of the magnitude Mmax.. The respective probability density function (PDF) of earth- quake magnitude for the exponential distribution is given as [Page, 1968]: fM (m) = k ? exp (-? (m-Mmin)) for Mmin ? m? Mmax (3.9a) where k= [1-exp (-? (Mmax-Mmin))]-1 fM(m) = 0 for m < Mmin and m > Mmax (3.9b) The characteristic parameters defined for each seismic source zone are: Average rate of occurrence or mean seismic activity rate ? (which is the parameter of the Poisson distribution), Level of completeness of the earthquake catalogue Mmin, 38 Maximum possible earthquake magnitude Mmax, Gutenberg-Richter parameter b (which indicates the relative number of large and small earthquakes, ?= b ln10) , Focal depth, and Regional attenuation relationship for the strong ground motion. 3.2.5.1 b-value and activity rate ? The b-value is expected to be regionally stable with variations less than the uncertainty limits, while the activity rate ? is liable to vary substantially from one seismic source zone to another. These parameters were obtained using the Seisan 8.1 software BVALUE program [Havskov and Ottem?ller, 2006] using a Nordic input file with magnitude increments ?m= 1.0, 0.5, 0.25 and 0.1, on the data selected from the earthquake catalogue corresponding to each seismic source zone. The program gives the number of events selected between magnitude Mmin and Mmax, cumulative number of events and duration of the catalogue in years. An example of earthquake magnitude recurrence relations for the source zone 1, 2 and 3, using magnitude increment equal 0.25, are shown in fig.3.5a-b-c. The total number of events of magnitude greater or equal to 4 used to calculate the b-value for each seismic source zone was 99 for Congo Basin, 99 for Upemba and Moero Rifts, and 624 for the Western Rift zone. The maximum standard deviation in the calculation of b-value in these three seismic source zones was ?0.1 39 Figure 3.5- Cumulative magnitude-frequency plot for a magnitude increment ?m= 0.25 for the three main seismic source zones used in seismic hazard analysis: (a) Upemba?Moero Rift, (b) Congo basin, and (c) Western Rift Valley. 40 3.2.5.2 Maximum magnitudes At present, there is no generally accepted method for estimating the value of Mmax. Mmax is sometimes inferred through other available information, such as geology, paleoseismicity or subjective judgement of the scientist. If the above information is not available, it is usually set to half a magnitude unit higher than the maximum observed magnitude. More physically-based estimates are made using functions that relate magnitude to rupture length (for faults), or to the length of tectonic features (for area sources). Also analogies with regions with similar tectonic history can give insight on the maximum magnitude [McGuire, 1993, Wang and Tim Law, 1994, Havskov and Ottem?ller, 2006]. As we don?t have well-documented information about faults or the slip and displacement of paleoseismic events, we used a procedure for the evaluation of Mmax developed by Kijko and Sellevoll [1989] and Kijko [2004] in this study. This procedure is free from subjective assumptions and is dependent only on seismic data. The procedure can be applied in the extreme case when no information about the nature of the earthquake magnitude distribution is available. In this procedure, it is assumed that: the value of the magnitude Mmin is known, the magnitudes of all the main shocks that occurred in the area within a specified time period T are independent and identically and randomly distributed with cumulative distribution function (CDF), FM(m). The maximum earthquake magnitude Mmax (an unknown parameter) is the upper limit of the range of magnitude. If Mmaxobs is the largest observed magnitude, the maximum regional magnitude Mmax can be approximated by the formula [Kijko and Graham, 1998] 41 dmmFMM n M M M obs )]([maxmax max min (3.10) For the Gutenberg-Richter frequency-magnitude relation, the respective CDF of magnitudes bounded above by Mmax is [Page, 1968]: 0)(mFM for m< Mmin min)]max(exp[1 min)](exp[1 )( MM M mFM for Mmin ? m ? Mmax (3.11) 1)(mFM for m> Mmax where ?= b ln(10), and b is the b-parameter of the Gutenberg-Richter relation. Following (3.10) and using the approximation of Cramer [1961], the estimator of Mmax for the frequency-magnitude distribution is obtained as a solution of the equation: )exp(min )exp( )(1)(1 maxmax 2 12 nM n nEnE MM obs (3.12) where min)]}max(exp[{ MM1 nn1 min)]max(exp[ MM nn 12 and E1 (.) denotes an exponential integral equation. The function E1 (.) is defined as z dzE /)exp()(1 , and can be approximated as 42 )exp( )( )(1 21 2 21 2 z bzbzz zazaz zE (3.13) Where a1= 2.334733, a2= 0.250621, b1=1.681534 and b2= 1.681534 [Abramowitz and Stegun, 1970]. It must be noted that in its current form, equation (3.12) does not constitute an estimator for Mmax since expressions n1 and n2, which appear on the right-hand side of the equation, contain also Mmax. The assessment of Mmax can be obtained by the iterative solution of equation (3.12). From equations (3.10) and (3.12), the approximate variance of the maximum regional magnitude Mmax is 2122 )]exp(min )exp( )(1)(1 [max)( nM n nEnE MVar M (3.14) where ?M 2 is the variance in the determination of the largest observed magnitude Mmaxobs. Equation (3.12) has been used to estimate the maximum possible earthquake magnitude in several seismically active areas such as China [Yurui and Tianzhong, 1997]. We developed a program using the iterative solution method and calculated Mmax by taking a starting value of Mmax as Mmax=Mmaxobs+0.5. The maximal standard deviation in the calculated Mmax using formula (3.14) was ? 0.5. 43 3.2.5.3 Focal depths Due to the poor coverage of seismographic stations in the region, the determination of focal depth is generally poor for events that occurred outside our seismographic network. However, micro-seismic studies indicate that, in general, the depth of earthquakes foci in the Western Rift Valley of Africa ranges between 10 and 20 km (Zana [1977], Zana and Hamaguchi [1978], De Bremaecker [1959] and Wohlenberg [1967, 1968]. In this study, an average value of 15 km with 5 km uncertainty was considered. 3.2.5.4 Attenuation relations As no strong ground motion database is available to use in a regression analysis to estimate ground motion parameters as a function of magnitude and distance, it was necessary to resort to other, less direct, data sources. Jonathan [1996] and Twesigomwe [1997] attempted to establish an average attenuation relation for the region. Jonathan?s relation is based on random vibration theory using some earthquakes recorded by the digital stations in the region. Twesigomwe?s relation is a modification of the previously established relation by Krinitzky et al. [1988] using regional shear-wave velocity and Q values determined by other workers like Gumper and Pomeroy [1970] in the Uganda Craton. The two relations are given below: ln a =3.024 + 1.030 Mw -1.351 ln R ? 0.0008 R ? ? [Johathan, 1996] (3.15) ln a= 2.832 + 0.866 Ms ? ln R ? 0.0025R ? ? [Twesigomwe, 1997] ; (3.16) where a is the ground acceleration (cm/sec2), 44 R is the hypocentral distance (km), and ? is the error term. The standard deviation of both these relations is ?=0.6. Both relations were developed for hard rock conditions. Stochastic methods has been widely used to predict site-specific time histories of strong motion for future large earthquakes (e.g. Boore [1983]; Hlatywayo [1995]), but they generally do not incorporate the response that is observed for a specific site, or the change in frequency content with time in seismograms caused by the late arrival of surface waves and attenuated shear waves. Mavonga [2007b] derived an attenuation relationship based on the simulation of the strong motion of large earthquakes using recordings of small earthquakes [Frankel, 1995; Irikura, 1983, 1986]. Recordings of small earthquakes adjacent to the expected large earthquakes have been treated as an empirical Green's function (EGF) between the source and that site [Hartzell, 1978]. These recordings contain the site response, scattered waves, surface waves and other path effects, and can be summed to produce realistic time histories at that site for a large earthquake. On the basis of the assumption of identical seismic response at three hard rock sites (Lwiro, Katale and Kunene) located within the Western Rift Valley of Africa and the relationship between the fault area and moment magnitude of the earthquakes derived from aftershock studies of major earthquakes that occurred in the Kivu Province from 1994 to 2002 [Mavonga, 2007a], an attenuation relationship was calculated based on simulating strong motions of large 45 earthquakes using recordings of small earthquakes (aftershocks). Results showed that the attenuation of mean peak ground acceleration (Y) was related to the moment magnitude (Mw) and the epicentral distance (R), according to the following power laws: Mw=5.0 Y= 1.42 exp (1.43 Mw) R-1.1 Mw=5.5-6.5 Y= 1.42 exp (1.43 Mw) R-1.2 (3.17) Mw=7.0 Y= 1.42 exp (1.43 Mw) R-1.3 Mw=7.5 Y= 1.42 exp (1.43 Mw) R-1.4 The peak ground acceleration (Y) and epicentral distance (R) are expressed in gal and kilometre (km), respectively. In this study, we considered the relations (3.17) under the following form: ln Y = c1 + c2 Mw + c3 ln (R + C4 ) + C5 R + ? ?? N(0, ? ) (3.18) R is the epicentral distance (in km) and the random error ? has a normal distribution with mean 0 and variance ?2. The parameters of the model are: c1= -6.53857, c2 =1.43, c3= -1.5, c4=0, c5=0 and ? =0.70. In equation (3.18), Y is expressed in unit of g. Note that ? is the uncertainty in the estimated maximum value of PGA. It includes the value of ?ln(PGA) due to the random scatter of the logarithm of acceleration, uncertainty arising from the erroneous determination of earthquake magnitude (?M 2) and distance (?R 2). It was determined by means of equation [Tinti and Mulargia, 1985]: ? = (?ln(PGA) 2 + c2 2 ?M 2 + ?R 2 ( c5 + c3/R) 2 )1/2 (3.19) 46 We used two epicentral distances (R=25 km and R=100 km) and assumed the maximum uncertainty values in the determination of magnitude and distance to be ?M=0.4 and ?R =10 km, respectively. The maximum value of ?ln(PGA) was 0.38. To select ground-motion estimation equations for input in the seismic hazard calculation, we considered five attenuation relations for hard rock conditions, of which three were derived from data for central and southern Africa: 1. Mavonga [2007b] from data recorded in the Western Rift Valley, 2. Twesigomwe [1997] for Uganda, 3. Jonathan [1996] for eastern and southern Africa, 4. Atkinson and Boore [2006] for eastern North America, and 5. Somerville [2001] for central North America. The North American attenuation equations are the best known and most frequently used attenuation formulas for stable continental regions (e.g. Kijko et al., 2002). Our choice in adopting these equations was motivated by the following considerations: Sub-Saharan Africa is mostly a stable intra-plate region characterized by a relatively low level of seismic activity, with earthquakes randomly distributed in space and time. The only parts of Sub-Saharan Africa that do not display the characteristics of an intra-plate region are the East African Rift System and the Cameroon volcanic line, where earthquakes are associated with active fault zones and volcanic activity (Kubanza et al., 2006). The Democratic Republic of Congo and surrounding areas encompass both cratonic regions and the rift zone. Based on our 47 knowledge of the investigated area, it is not yet possible to delimit exactly the boundary between the craton and rift zone due to the sparse distribution of seismographic stations in the region and poor knowledge of the deep earth structure of the study area. The dependence of ground motions on geographical region is an important consideration when selecting an attenuation model. However, these questions can be particularly difficult to tackle in regions of the world where little observed strong-motion data is available, since there are few records to validate the choice of the model. Thus it is often preferable to use well-constrained models, based on data from a other similar regions, than to predict motions using local models that are poorly-constrained (Douglas, 2007). It is recommended that hazard be represented by several curves with the median or mean curve rather than a single curve. In this way both the uncertainty and the central value of the hazard are represented and may be considered for mitigation decision (Risk Engineering, Inc., 2008). The five attenuation curves are shown in Figure 3.6. The figure shows a close agreement between the relations of Jonathan [1996], Mavonga [2007b] and Atkinson-Boore [2006]. The relation of Twesigomwe [1997] is close to that of Somerville [2001], but both present a large departure from the first three. Therefore, in the computation of seismic hazard, we used the Jonathan [1996], Mavonga [2007b] and Atkinson-Boore [2006] attenuation relationships. The Mavonga [2007b] attenuation equation, obtained from the Western Rift Valley data, attenuates more strongly than those derived from cratonic or stable region 48 data. This is attributed to the presence of heterogeneous material generated either by the volcanism or rifting in the shallow part of the crust. Figure 3.6 Attenuation curves obtained for Central and Southern Africa data, together with the relations for Eastern North America [Atkinson and Boore, 2006] and Central USA [Somerville et al., 2001]. These curves are for magnitude Mw7.5. The PGA in the Atkinson and Boore [2006] and Somerville et al. [2001] relations was obtained at a fixed frequency of 50 Hz in the spectral acceleration spectra curve. 49 3.2.6 Input parameters for hazard computation The calculated parameters used in the seismic hazard computation are listed in Table 3.1. Table 3.1: Seismicity parameters for the area source zones SOURCE ZONES Mmin Mmax Beta(?) Lambda(?) 1.Upemba ?Moero Rift 4 6.99 2.07 1.95 2. Congo Basin 4 6.09 2.39 1.81 3.Western Rift 4 7.79 1.84 6.46 3a. Ruwenzori area-Lake. Edouard trough 4 7.29 1.84 2 3b. Virunga volcanic complex-Rutsuru-Masisi 4 5.99 1.84 1.13 3c. Lake Kivu basin ?Ngweshe-Ruzizi 4 6.67 1.84 0.74 3d. Tanganyika Rift 4 7.79 1.84 3.17 Mmin: lower bound magnitude; Mmax: maximum expected upper bound magnitude; b-value: slope of magnitude-frequency relation; Beta (?): ln(10)x b-value, Lambda (?): annual number of earthquakes above the lower bound magnitude 50 3.2.7 Computation of seismic hazard 3.2.7.1 Theoretical background The seismic hazard calculation at a given site due to multiple earthquake source zones can be represented by the following equation [Reuter, 1990] which uses the total probability theorem to calculate the probability of a ground motion (e.g. Peak ground acceleration, velocity or displacement) being exceeded at a given site. n i iRMi dmdrmrMfmfrmaAPaH ii 1 ),(|)(],|[)( (3.20) where H (a) is the annual frequency of earthquake that produces ground motion with amplitude A greater than a. is the annual rate of earthquakes (with magnitude greater than Mmin) in the i-th source area. P [A>a | m, r] is the probability that an earthquake of magnitude m at distance r produces a ground motion amplitude A at the site that is greater than a. It is obtained from the cumulative lognormal distribution with specified standard deviation ? as: duyurmaAP a ] 2 1 )(exp[ 1 2 1 1],|[ 2 ln 2 (3.21) Where y = C1 + C2 M + C3 ln (R + C5) + C4 R +? ?? N (0, ?) for area sources. C1 , C2. C3, C4 and C5 are empirical constants M is the magnitude of earthquakes 51 ? is a random error which has a normal distribution with mean 0 and variance ?2. )(mf iM is the probability density function of earthquake magnitude. It depends upon the earthquake magnitude recurrence model. Often (m) is assumed to be an exponential truncated at lower and upper limits by Mmin and Mmax as expressed in equation (3.9). ),(| mrMf iRi is the probability density function of earthquake-site distance. It depend upon the geometry of earthquake source, which usually takes the form of a point, line or bounded surface. The occurrence of earthquake is assumed to follow the Poisson probability density function. Therefore, the probability of exceedance r(a) of the ground motion a is often expressed by the following equation: ))(exp(1)( aTHar (3.22) Where T is the time period (number of year) for which we want to know the probability; H(a) is the annual rate of exceedance of ground motion a, and 1/H(a) is the return period. 3.2.7.2 Seismic hazard analysis The probabilistic seismic hazard for DRC and surrounding areas was computed using the software EZ-Frisk package (Version 7.24) [Risk Engineering Inc, 2007]. The inputs to probability hazard assessment are generally uncertain because they are based on the interpretation of limited data. This study on the DRC and 52 surrounding areas is no exception. To account for some of these uncertainties, modern seismic hazard analysis techniques regard input parameters as random variables. Their uncertainties are accounted for using a logic tree approach, which accommodates alternative input parameter values for a range of hypotheses, which are individually weighted and whose contributions to seismic hazard are separately evaluated and statistically combined (Mc Guire, 1993). However, EZ-Frisk does not use a logic tree directly. To treat epistemic uncertainty in the input to the analysis, it allows the use of multiple earthquake magnitude recurrence models (see glossary) after defining multiple identical area sources with probability of activity less than one in each seismic source zone. The EZ-Frisk program can also calculate seismic hazard using multiple attenuation equations and depths. EZ-Frisk presents results as a total hazard curve, uniform hazard spectra, contribution of hazard by seismic source and deterministic spectra. The total hazard curve shows the annual frequency of exceedance as a function of spectral acceleration for each attenuation equation for a particular spectral period. The uniform hazard spectra show the spectral acceleration as a function of spectral period for a specified return period. The deterministic spectra show spectral acceleration as a function of spectral period for a specified fractile of the attenuation dispersion. In this study, EZ-Frisk was used to compute the mean peak ground acceleration for two combinations of seismic source zones, three attenuation equations and two alternative focal depths at the site of interest 53 3.3 Results 3.3.1 Seismic hazard map The input parameters used for seismic hazard computation are listed in Table 3.1. The output is a statistical estimate of the annual chance of exceedance as a function of PGA. A 0.5 degree grid was used to determine the PGA at various sites for a specific return period. The values obtained were used to create the maps shown in Figures 3.7 (a), (b) and (c). These figures show the mean PGA value (in unit of g) for a 2%, 5% and 10% probability of exceedance in 50 years which correspond to return period of 2475, 975 and 475 years, respectively. The highest PGA values are found in the Lake Tanganyika Rift zone where the PGA values of 0.32g, 0.22g and 0.16 g are expected to be exceeded with probability 2%, 5% and 10% in 50 years, respectively. The regions with the next high level of hazard are the other sub-zones in the Western Rift Valley, with the exception of the Virunga volcanic complex-Rutsuru-Masisi sub-zone where the seismic hazard (due mainly to the volcanic activity) is only moderate. The seismic hazard in the Congo basin diminishes with distance from the Western Rift Valley until, at a distance of about 450 km, the chance of exceeding a PGA of 0.05g (the threshold value of engineering interest) is less than 10% in 50 years. 54 Figure 3.7 (a): Distribution of mean PGA values (in unit g) in the DR Congo and surrounding areas computed for 2 % chance of exceedance in 50 years. The filled circles and triangles indicate the site PGA values and main cities, respectively. KIS=Kisangani, BNA=Bunia, BTC=Butembo, GOM=Goma, BUK=Bukavu, KAL=Kalemi, LUB=Lubumbashi, MBJ=Mbuji-Mayi, KND=Kindu, and KNG=Kananga. The geographic coordinates of the cities are listed in Table 3.2. 55 Figure 3.7 (b): Distribution of mean PGA values (in unit g) in the DR Congo and surrounding areas computed for 5 % chance of exceedance in 50 years. 56 Figure 3.7(c): Distribution of mean PGA values (in unit g) in the DR Congo and surrounding areas computed for 10 % chance of exceedance in 50 years. 57 3.3.2 Total hazard curve at the sites of the cities Figures 3.8 (a)-(n) show the three total hazard curves (one for each attenuation equation) and the mean total hazard curve at each of the following sites (in the DRC unless specified otherwise): Bujumbura (Burundi), Bukavu, Bunia, Butembo, Goma, Kalemi, Kananga, Kigali (Rwanda), Kigoma (Tanzania), Kindu, Kisangani, Lubumbashi, Mbuji-Mayi and Uvira. These sites are shown in Figure 3.7 by filled triangles. The PGA exceeding indicated values with probability of 2%, 5% and 10% in 50 years was estimated for these sites using both a probabilistic and a deterministic approach. These values are listed in Table 3.2. The probabilistic approach is dependent on the location of the site relative to all seismic source zones, while a deterministic approach considers only the largest possible magnitude event that is expected to occur in the seismic source zones closest to the site location. The following levels of hazard were found: Very high in the Tanganyika Rift region, where the cities of Bujumbura (Burundi), Uvira, Kalemi and Kigoma (Tanzania) are located. The PGA, averaged for the four cities, to be exceeded with probability of 2%, 5% and 10% is 0.32g, 0.20g and 0.15g, respectively. High in the Lake Kivu Basin and Ruwenzori areas, where the cities of Bukavu, Goma, Butembo, Bunia and Kigali are located. An average PGA in excess of 0.27g, 0.18g and 0.13 g with probability of 2%, 5% and 10%, respectively, was determined in these areas. Moderate in the Upemba and Lake Moero area and in the part of Congo basin close to the Western Rift. This area includes the cities of 58 Lubumbashi, Mbuji-Mayi and Kindu, An average PGA in excess of 0.13 g, 0.08g and 0.05 g with probability of 2%, 5% and 10%, respectively, was determined in these areas. Low at Kananga and Kisangani. An average PGA in excess of 0.07g, 0.05g and 0.03 g with probability of 2%, 5% and 10%, respectively, was determined in these areas. Table 3.2: Peak ground acceleration at selected cities determined using probabilistic and deterministic approaches for two seismic source zone models (S1 and S2) of the Western Rift Valley CITY LATITUDE LONGITUDE 2% 5% 10% S2 S1 S2 S1 S2 S1 S2 S1 BUJUMBURA -3.378 29.363 0.3155 0.319 0.1947 0.1989 0.1311 0.135 2.416 2.416 BUKAVU -2.53 28.85 0.274 0.3185 0.1777 0.198 0.1239 0.1337 0.771 2,416 BUNIA 1.528 30.255 0.2728 0.3189 0.1708 0.1989 0.1169 0.1349 1.435 2.416 BUTEMBO 0.138 29.288 0.2724 0.3187 0.1699 0.1983 0.1157 0.1342 1.435 2.416 GOMA -1.683 29.231 0.2709 0.319 0.1782 0.1988 0.1266 0.1347 0.771 2.416 KALEMI -5.93 29.176 0.3196 0.3187 0.199 0.1984 0.1346 0.1342 2.416 2.416 KANANGA -4.083 21.7 0.07567 0.09048 0.0455 0.0522 0.02993 0.03342 0.4315 0.4315 KIGALI -1.953 30.059 0.2705 0.3191 0.1744 0.1991 0.1215 0.1352 0.771 2.416 KIGOMA -4.796 30.402 0.3197 0.319 0.1994 0.199 0.1353 0.1351 2.416 2.416 KINDU -2.921 25.915 0.07719 0.09239 0.04815 0.05579 0.03322 0.03793 0.4315 0.4315 KISANGANI 0.425 24.029 0.0757 0.09059 0.04566 0.05259 0.03025 0.03424 0.4315 0.4315 LUBUMBASHI -11.669 27.485 0.1259 0.1259 0.07776 0.07777 0.05227 0.05229 1.06 1.06 MBUJI MAYI -6.16 22.885 0.1186 0.1201 0.0712 0.07217 0.04607 0.0468 1.06 1.06 PROBABILISTIC APPROACH DETERMINISTIC APPROACH PGA(g)PGA(g) T1=2475 YEARS T2=975 YEARS T3=475 YEARS 59 (a) (b) Figure 3.8 Total hazard curve at the site (a) Bujumbura (b) Bukavu 60 (c) (d) Figure 3.8 Total hazard curve at the site (c) Bunia (d) Butembo 61 (e) (f) Figure 3.8 Total hazard curve at the site (e) Goma (f) Kalemi 62 (g) (h) Figure 3.8 Total hazard curve at the site (g) Kananga (h) Kigali 63 (i) Figure 3.8 Total hazard curve at the site (i) Kigoma (j) Kindu (j) 64 (k) (l) Figure 3.8 Total hazard curve at the site (k) Kisangani (l) Lubumbashi 65 (m) (n) Figure 3.8 Total hazard curve at the site (m) Mbuji Mayi (n) Uvira 66 Figure 3.9 shows the contribution of each seismic source zone in the study area to the total hazard at Bukavu site. The total hazard at this city is mostly dominated by the hazard due to seismic source zone 3c (Lake Kivu Basin-Ngweshe-Ruzizi) because Bukavu city lies within this seismic source zone. By using eq. (3.22), this chart shows the PGA of 0.05 g (the threshold value of engineering interest), corresponding to an annual frequency of exceedance of 0.01, is expected to be exceeded with probability of 10%, 39% and 63% in 10 years, 50 years and 100 years, respectively. Figure 3.9 Contribution of hazard by seismic source at Bukavu site 67 3.4 Discussion 3.4.1 General overview The main steps and parameters of the probabilistic seismic hazard assessment conducted in this study can be summarized as follows: 1. Build a catalogue from instrumental seismic data with a unified magnitude that provides information on the location and frequency of earthquake occurrence during the past years. 2. Delineate seismic source zones based on geological and seismological evidence. These source zones describe the potential locations of future earthquakes within the study area. 3. Evaluate, for each seismic source zone, earthquake seismic parameters such as maximum expected magnitude, activity rate and b-value of the Gutenberg-Richter relation. 4. Predict future ground motion using an appropriate regional attenuation relationship for the strong motion for the study area between magnitude, distance and site conditions. 5. Compute seismic hazard (the probability that a specified ground motion level at a site will be exceeded during a particular time period), using the above parameters as input, to characterise each seismic source zone. 68 3.4.2 Source of uncertainties The input parameters used in seismic hazard computation are subject to several uncertainties: 1. The results obtained in this study are based on a period of about 90 years duration, and assumed to be complete for events of magnitude Mw>4. However, African plate boundaries are generally characterized by slow relative motions (?2mm/yr to ?15mm/yr). Hence, large earthquakes have extremely long recurrence times [Hartnady and Benouar, 2007; DeMets et al., 1990; Walpersdorf et al., 1999; Lowrie, 1997; Asfaw et al., 1992, Jestin et al., 1994, Chu and Gordon, 1999, Bilham et al., 1999]. As the maximum credible magnitude was estimated on the basis of the observed maximum magnitude, this important parameter could be underestimated. 2. The regression relations used to homogenise magnitude introduces errors that are later incorporated in the calculation of the b-value and the activity rate. 3. Observations of regional strong ground motion are lacking. The attenuation relations used in the analysis are based on numerical simulations. 4. The delimitation of seismic source zones also requires significant interpretation on the part of the analyst as the seismic stations are sparsely distributed and the geology is poorly documented. 69 3.4.3 Assessment of the sensitivity of the outputs to input data To assess the sensitivity of the outputs to these uncertainties, we considered two alternative source delimitations, three attenuations equations and two different source depth estimates (10 km and 20 km). To determine the effect of delimitation of seismic source zones on estimates of the PGA at selected sites, both probabilistic and deterministic approaches were used, considering two delimitations of the Western Rift Valley (WRA) seismic source zone. In the first scenario (denoted S1), the Western Rift Valley seismic source zone was considered to be one structure with a single rate of activity, b-value and Mmax. In this model, every point in the Western Rift Valley is a candidate for the largest magnitude event. In the second scenario (denoted S2), the Western Rift Valley seismic source zone is considered to be segmented into 4 sub-zones. The rate of activity ? and the maximum magnitude Mmax may differ from one sub-zone to another, but the b-value is considered to be the same for all these sub-zones as they are all governed by the same tectonic setting. The results are listed in Table 3.2 Using the probabilistic approach, the PGAs calculated using the S1 model were slightly greater (by about 0.01g) than those produced by the S2 model. This is not significant from the engineering point of view. Therefore, in this study, we adopted the S2 model in the final computation of seismic hazard because, on the basis of the seismicity and the geology of the region, we believe that the S2 model is more realistic. The probability of observing an earthquake with magnitude greater than 7 in some sub-zones (e.g. the Virunga volcanic complex and vicinity) 70 is small due to the high heat in the crust generated by volcanic activity. It may be impossible to accumulate sufficient strain to generate this maximum magnitude. In the case of deterministic approach, the maximum difference in the strong ground motions predicted by the two scenarios was more than 1g (Compare Bukavu, Goma and Kigali, Table 3.2). Using S2 model and deterministic approach, the maximum credible PGA to be exceeded determined for the site Bukavu, Kalemi, Butembo and Lubumbashi are 0.771g, 2.416g, 1.435g, 1.06g, respectively. However, the probabilistic hazard curves of the respective sites (Figures 3.8(b), (f), (d) and (l)) show that the annual probability of such events occurring close enough to Bukavu, Kalemi, Butembo and Lubumbashi are 0.00002, 0.000004, 0.000004 and 0.000002, respectively. These equate to events with 50 000, 250 000, 250 000, 500 000 year return period, respectively. Thus, the probabilities that such ?scenario earthquakes? will occur are extremely small, but not null. The comparison of PGA values obtained using these two approaches show large differences. The deterministic approach is essentially based on a worst case scenario where the event occurs close to the site. However, the approach does not provide any indication of how likely this event is to occur during the lifetime of the structure. The deterministic approach is usually used for significant structures like power plants, large dams, large bridges, etc. In developing such scenario earthquakes, considerations should be given to the potential fault locations which can generate these earthquakes. These two methods can complement one another to provide additional insights to the seismic hazard or risk problem. One method 71 will have priority over the other, depending on the seismic environment, the scope of the project and how quantitative are the decisions to be made. As the locations of faults that may potentially cause the scenario earthquake are not well documented, we opted to produce a hazard map based only on the probabilistic approach. This quantification of the seismic hazard is based mainly on the seismic and geological history available in the area. Incorporation of any new specific information may yield a different result. Source zonations based on new paleoseismic information may improve the detailed hazard map. Also, site response studies will provide improved local hazard curves. 3.4.4 Comparison with previous work in the Western Rift Valley The results obtained in this study were compared with previous work in the Western Rift Valley of Africa. Midzi et al. [1999] presented a seismic hazard map for the Eastern and Southern Africa. Their study covered a large area of approximately 50? x 50?. At such a large scale, only regional structures were considered when delimiting source zones. They reported a high PGA in excess of 0.22 g (240 gals) in Southern Sudan, Western Rift and Northern Tanzania with probability of 10% in 100 years. Twesigomwe [1997] presented a probabilistic seismic hazard map of Uganda. He reported a PGA in excess of 200 gals (0.20g) with probability of 5% in 50 years in or close to the Western Rift south of latitude 0.5?N. 72 Ferdinand [2007] analyzed the probabilistic seismic hazard in the city of Arusha in Tanzania, located about 200 km from the Lake Tanganyika area. This city is subject to threats from both the Western and Eastern Rift. He found exceedance of the PGA values of 0.15, 0.20 and 0.27g for return periods of 475, 975 and 2475 years, respectively. Using the parametric-historic procedure of Kijko and Graham [1998, 1999], Kijko et al. [2003] compiled a probabilistic seismic hazard map for Sub-Saharan Africa. Their map depicts a 10% probability of exceeding the calculated PGA at least once in 50 years. They found high hazard regions (areas that have a 10% probability of experiencing a PGA of 0.15 g and greater) in the south-western Cape region (South Africa), north-western tip of the Namibia-Angola border, parts of Mozambique, the Zimbabwe- Zambia border, and much of the East Africa Rift Valley. Our results are in general accord with all these earlier studies. 73 CHAPTER 4 VOLCANOGENIC SEISMICITY IN THE VIRUNGA AREA 4.1 Introduction The Virunga volcano group (Figure 4.1), which is located at the northern edge of Lake Kivu, consists of eight major volcanoes divided into three subgroups: the Eastern, the Central and the Western subgroups (Table 4.1). The volcanoes in the Western subgroup (Nyiragongo and Nyamuragira) have been the most active since 1882 (Demant et al., 1994, Clay and Cassandra, 1995). The Nyamuragira volcano (1.42?S, 29.2?E) is located 15 km northwest of Nyiragongo. It has exhibited two types of eruptive activity during the last century: lava lake activity in the summit caldera from 1921 to 1938; and flank eruptions that issued lava through new fissures, which usually began as lava fountains that gradually converged to form cinder cones. Since the Tsambene eruption (from 1938 to 1940), all eruptions until November 2006 have been flank eruptions. The lava lake has disappeared completely. The Nyiragongo stratovolcano (1.52?S, 29.2?E), about 20 km north of Lake Kivu, has been active since 1884 (Clay and Cassandra, 1995). Its eruptive activity is characterized by the appearance of a lava lake in the summit crater due to intermittent short pulses. It had a lava lake in the summit crater from 1928 to 1977. On January 10, 1977, approximately 22x106 m3 of extremely fluid lava erupted from flank fissures and flowed down slope at up to 60 km/hour, killing about 70 people and reaching within 600 m of the Goma airport. 74 Figure 4.1 Map showing the eight volcanoes of the Virunga volcanic zone. NY=Nyiragongo; NM=Nyamuragira; KA=Karisimbi; MI=Mikeno; VI=Visoke; SA=Sabinyo; GA=Gahinga; MU=Muhavura. The filled small circles and triangles indicates cities and seismographic stations, respectively 75 Table 4.1: Eight major volcanoes of the Virunga area Location Name of the volcano Altitude (m a. s. l) Eastern Group Muhavura 4127 Gahinga 3474 Sabinyo 3647 Central Group Visoke 3911 Karisimbi 4506 Mikeno 4437 Western Group Nyamuragira 3056 Nyiragongo 3470 On 17 January 2002, the Nyiragongo volcano erupted again, extruding about 20x106 m3 of lava. Lava flows originated out of the central crater as well as several locations along a huge system of fractures that rapidly developed along the entire southern flank of the volcano towards the city of Goma. Two main flows entered the town producing major devastation, and forced the rapid exodus of 300,000 to 400,000 people. About 15 % of the surface of the town was affected, including parts of the airport, most of the business centre, and the dwellings of 76 about 120,000 people. It was estimated that 150 people died as a consequence of eruption. This eruption was exceptional. It lasted only two days but it was accompanied by a swarm of more than 100 felt earthquakes during about 5 days. Several of these earthquakes were recorded regionally and teleseismically. The largest earthquake magnitude in this swarm was Mb5.2 [Sadaka et al., 2003]. Earthquake hypocenter determination is a powerful tool that may be used to assess magmatic activity and issue early warnings in the case of a prominent eruption. However, the problems of earthquake location and of velocity structure determination are interdependent [Chiarabba et al., 2000; Zhang and Thurber, 2003; Stephen, 2005]. In the computation of earthquake hypocenters in the Virunga area, it is desirable to use a P- and S-wave velocity structure specific for the Virunga area. However, knowledge of the structure of the crust and upper mantle in the Virunga area is poor compared with many other active volcanic regions of the world, particularly in Europe, Japan, North and South America. Apart from seismic reflection studies of the shallow structure (0-10 km) carried out in the Lakes Tanganyika, Malawi and Rukwa by the PROBE project [Rosendahl et al., 1992], data on deeper structure in the Western Rift are extremely sparse [e.g. Prodehl and Mechie, 1991]. No deep seismic reflection or refraction data is available for the Western Rift near the Virunga volcanic area. Nolet and Mueller [1982] used teleseismic data to derive a regional model of the Western Rift lithosphere. Their model has to be considered as a representative of a 77 broad area (1?S to 10?N latitudes) and contains a 35 km thick crust and relatively low upper mantle velocities to a depth of 140 km. In this study, receiver function analysis of teleseismic events recorded by two broadband stations were used to study the crustal structure beneath the Virunga area. The new velocity model was then used to determine the hypocenters of earthquakes that occurred in the period 18 August 2002 to 7 May 2004 prior to 8 May 2004 eruption of Nyamuragira, and 1 July 2004 to 27 November 2006 prior to the 27 November 2006 Nyamuragira eruption. 4.2 Data acquisition The seismograms used in this study were provided by the seismological network of Goma Volcano Observatory (GVO). The observation network is composed of seven stations: Katale (KTL), Luboga (LBG), Kunene (KNN), Rusayo (RSY), Kibumba (KBB), Goma and Bulengo (BLG). These stations are each equipped with a short-period Kinemetrics vertical-component SS-1 ranger seismometer (To=1 sec) connected to a PS-2 portable seismic recorder instrument. Signals from the seismometer are amplified and filtered in the amplifier module. The overall maximum of the PS-2 is 1 mm deflection for 1?V of input voltage. The amplifier panel has controls for amplifier gain and filter setting. The low-pass filter has -12 dB/ Oct roll-off and can be set to 2.5, 5.0, 12.5, 25.0 and 50.0 or out (70 Hz). The high-pass filter with +12 dB/ Oct roll-off can be set at 0.1 Hz or 5 Hz. At our stations, the low-pass and high-pass filter were set to 12.5 and 0.1 Hz, respectively, and the amplifier gain varied from 36 to 66 dB according to the response of the site. 78 On November 2003, GVO began the deployment of three-component short-period Lennartz (LE-3D/5sec) seismometers at these stations. On May 2004, the Lennartz short-period seismometer at KNN and KBB was replaced with Nanometrics Trillium 40 broadband seismometers. Signals from these stations are locally digitized from a modular acquisition system (GAIA), specifically developed by the ?Istituto Nazionale di Geofisica e Vulcanologia? (INGV), with a sampling frequency of 50 Hz and an A/D resolution of 24 bits, and are telemetered to the Goma base station where they are recorded in triggered and continuous files. 4.3 Crustal structure beneath two seismic broadband stations in the Virunga Volcanic area 4.3.1 Method 4.3.1.1 Receiver function By inverting receiver functions, it is possible to constrain the Moho depth and average shear wave velocity beneath a recording station [Owens et al., 1987; Langston, 1989; Ammon, 1991; Ammon and Zandt, 1993; Mangino et al., 1993; Gurrola et al., 1994; Kaspar and Ritter, 1998; Sandvol et al., 1988; Midzi and Ottem?ller, 2001]. 79 In the receiver function technique, the teleseismic body waveforms are used to image the crustal structure underneath isolated seismic stations. These waveforms contain information related to the source time-function, propagation effect through the mantle, and local structure underneath the recording site. The resulting receiver function is obtained by removing the effects of the source and mantle path. The receiver function method relies on the conversion of the P-wave signals from distant events incident to a discontinuity in the crust or upper mantle to an S- wave, which arrives at the station within the P-wave coda directly after the direct P-wave. As the S-waves travel more slow than the P-waves, a direct measure of the depth of the discontinuity is calculated using the difference in the arrival times of the direct P-wave and the converted phase (Ps), provided the velocity model is known. The phase obtained by Ps conversions, which is sensitive to the shear velocity contrast, has much stronger amplitude on the horizontal component than on the vertical component. In addition to the direct Ps converted phase, the multiple reflection phases resulting from the reverberation between the discontinuity and free surface are also seen on the receiver function traces. This method was presented in earlier studies by Langston [1979], Ammon et al. [1990] and Cassidy [1992]. 80 In this study, we used the source equalization procedure developed by Langston [1979] to remove the effects of near-source structure and source time-function. In the frequency-domain deconvolution, we used a water-level stabilization method and a low-pass Gaussian filter to remove any high frequency noise not filtered by the water-level method. The radial component of the receiver function H (?) is obtained by deconvolving the vertical component Z(?) of an earthquake recording from the radial component R (?), i.e. )()( )()( )( * G ZR H (4.1) Where ?( ? ) = max { Z ( ? ) Z* (? ), c max [ Z ( ? ) Z* (? ) ] } (4.2) Z* (?) is the complex conjugate and c is the water-level constant. The water-level constant is the limit value of the denominator in Equ. (4.1) to avoid dividing by 0 [Clayton and Wiggins, 1976; Mangino et al., 1993]. Typical values used are 0.0001, 0.001, 0.01 and 0.1. The water-level constant is chosen by examining the results of several trials and choosing the lowest water-level that produces acceptable noise levels in the corresponding receiver function. In this study, we used the value of 0.1. G(?) is the low-pass Gaussian filter commonly used to reduce high frequency noise in the receiver function. It is expressed by: } 4 exp{)( 2 2 a G (4.3) Where a is the Gaussian width factor. In this study, we used a=5.0 81 The tangential component for the receiver function is reproduced by the same procedure. For a laterally homogeneous structure, it should theoretically be zero. This provides a way to judge the degree of lateral variation in earth structure beneath a station (e.g. the presence of a dipping velocity interface). In receiver function studies with epicentral distances greater than 30?, it is acceptable to assume that the incoming P-wave is a plane wave. Assuming that the ray parameter p for the direct P-wave and Ps conversions is the same and H the thickness of the crust, the delay time of the Ps converted wave can be calculated [Kind and Vinnik, 1988] as follows: ))()(( 2 1 222 1 22 pVpVHTT pppPs (4.4) Similarly, the travel time delay for crustal reverberations such as PpPs (two P legs and one S leg) can be given as: ))()(( 2 1 222 1 22 pVpVHTT pspPpPs (4.5) The receiver functions are traditionally inverted to produce a model of the S-wave velocity structure under a given seismic station. The method used in this study was developed by Owens [1984]. It incorporates minimum roughness constraints to remove the non-uniqueness problem from each individual inversion. A singular-value decomposition method was used to compute the matrix inverse and solve the inversion problem. There is no guarantee that a unique inversion result will be obtained, as the method seeks to minimize the differences between observed and synthetic receiver functions. Therefore, using a priori information available for the area under investigation is helpful in constraining the solutions. 82 An important assumption of the time domain inversion technique is that the initial model must be close to the true earth velocity structure [Ammon et al., 1990]. The inversion procedure consists of preparing the observations, constructing an initial model, choosing a smoothness weight parameter, inverting the waveforms and assessing the significance of features in the results. The inversion is a search for models fitting the observations using a gradient-based algorithm. The accepted solution is that which corresponds to the best fit between observed and synthetic receiver functions. In this study, we used three initial velocity models: i. Nolet and Mueller [1982] model, ii. the average model of Bram [1975], Bonjer et al. [1970], and Nolet and Mueller [1982] models, and iii. a combination of (i) and (ii). The maximum number of iterations was set to 4, the maximum cubic perturbation to 0.75, the stopping perturbing velocity at 7.8 km/s, the maximum random perturbation equal 20% of cubic perturbation, the smoothness trade-off parameter to 0.1, the singular value truncation fraction to 0.01, the average horizontal slowness to 0.06 s/km and the Gaussian width factor to 5. The number of iteration required and the smoothest model were determined using the pre-processing computer programs ?snglinv? and ?smthinv?, respectively, provided in the package of computer programs ?RftnCodes? by Ammon [1997]. 4.3.1.2. Data Processing More than 100 teleseismic events were recorded at broadband station KNN and KBB for the period from May 2004 to December 2007. We chose 46 earthquakes with magnitude greater than 5.5 that have a good signal-to-noise ratio and located 83 between 30? to 92? from the stations KNN and KBB. The locations of hypocenter (latitude, longitude and depth) and magnitudes of these teleseisms were obtained from the USGS (Table A2). The located earthquakes and the stations KNN and KBB are shown in Figure 4.2. Figure 4.2 Epicenter distribution of earthquakes used to determine P receiver functions. The events have magnitude larger than 5.5 (mb) with epicentral distances between 30? to 92?. The broadband stations KNN and KBB are marked by the filled triangle. The filled lozenges, circles, squares and stars indicate earthquakes of magnitude range 5.5 to 6.0, 6.0 to 6.4, 6.4 to 6.6 and 6.6 to 8.6, respectively. Data was processed using the package of computer programs ?RftnCodes? provided by Ammon [1997]. As our data was originally sampled at 50 Hz, before deconvolution, we decimated it from 50 Hz to 12.5 Hz using the SAC command ?decimate with anti aliasing filter on?. We computed receiver functions using software ?pwaveqn? based on the source equalization procedure developed by Langston [1979]. Receiver functions for events with similar shape, ray parameter 84 and back azimuth were grouped at back azimuthal intervals of ?15? at each station [Midzi and Ottem?ller, 2001]. The events were mainly clustered in the northeast (NE) and southeast (SE) quadrants with a few in the southwest and northwest. In the NE quadrant, important clusters were observed between 30? and 45? and 75? and 90?. In the SE quadrant, a cluster was observed in the interval between 90? and 105?. A small cluster was observed between 270? and 285? in the NW quadrant. Some isolated events were also observed in the NE (interval 0-15?, 15- 30? and 60-75?), SE (interval 135-150?), SW (intervals 195-210?, 210-225?, 225- 240? and 240-255?) and NW quadrant (intervals 285-300?). Fig.4.3A-B (Left) illustrates stacked receiver functions of waveforms in the NE quadrant in the back azimuthal interval from 0-45? and 60-90? at the KNN and KBB stations. No events were observed between 45-60?. At KNN, in the azimuthal interval from 0-45?, two weak phases were observed at 2.75 sec and 5.25 sec after the direct P phase. At KBB, two weak phases were observed at 3.5 and 5.6 sec after the direct P phase. These two weak phases observed between 45- 60? may correspond to the P to S converted phases at a mid-crustal boundary (Conrad discontinuity) and the crust mantle-boundary (Moho). However, only one sharp phase was observed for the back azimuthal interval 60-90?, at 4.75 sec and 3.70 sec for KNN and KBB, respectively. The difference between the 0-45? and 60-90? receiver functions imply a significant variation in crustal structure with back azimuth. 85 As we are interested, at this first stage, on the first-order approximation of the seismic velocity structure in the Virunga area; we ignored back azimuthal variation and stacked all the receiver functions in the back azimuthal interval 0- 90? and obtained average receiver functions for the NE quadrant at KNN and KBB, which are shown in Figure 4.3A-B (Right), respectively. In the SE quadrant, at back azimuthal interval of 90? to 105?, two receiver functions provided by two deep earthquakes which occurred in Indonesia (Day 25/07/2004, O.T=14:35:19.0, Lat. =2.427 S, Long. =103.981, Mw=7.3, h=582 km and Day 08/08/2007, O.T=17:05:04.9, Lat. =5.859 S, Long. = 107.419 E, Mw=7.5, h=280 km) were analysed. Although their source depths differ significantly from other earthquakes in the group, these earthquakes are within the azimuth criteria of 15? and similar horizontal slowness (p is about 0.05 s/km). Comparison of their receiver functions with those of shallow events in the group does not show a large difference (Figure 4.4A-B). Amplitudes differ at KBB, though timing is similar. All the receiver functions in the SE quadrant, from back azimuthal interval between 90-105?, along with the isolated at interval 135-150?, were stacked to get an average receiver function at KNN and KBB representative of the SE quadrant. In the SW quadrant, only isolated events were observed at back azimuthal intervals 195-210?, 210-225?, 225-240? and 240-255?. These receiver functions were stacked. In the NW quadrant, a small cluster was observed in the back 86 azimuthal interval 270-285? and isolated events in the interval 285-300?. All these receiver functions were stacked. These average receiver functions representative of NE, SE, SW and NW quadrants at KNN and KBB are shown in Fig.4.5A and Fig.4.5B, respectively. 87 Figure 4.3 Stacked receiver functions of waveforms obtained at (A) KNN and (B) KBB stations in the NE quadrant at back azimuthal interval between 0? to 45?, 60? to 90? (Left) and 0? to 90? (Right) 88 Figure 4.4 Comparison of receiver functions of deep earthquakes and shallow events at back azimuthal interval 90? to 105? in the SE quadrant at (A) KNN and (B) KBB stations 89 Figure 4.5(A) Stacked receiver functions representative of NE, SE, SW and NW quadrants along with the average receiver function (AVG) for all quadrants computed at KNN station 90 Figure 4.5 (B) Stacked receiver functions representative of NE, SE, SW and NW quadrants along with the average receiver function (AVG) for all quadrants computed at KBB station 91 4.3.2. Results 4.3.2.1 Ps converted phases and Moho depth Based on results of previous crustal structure studies in the region (e.g. Prodehl et al., 1994; Dugda et al., 2005), we assume an average velocity of the crust in our study area to be =6.2 km/s and =3.6 km/s. Using the approximation of Kind and Vinnik [1988], we estimated the depth equivalent to the delay time between the direct P phase and converted phase Ps at KNN and KBB. The Ps-P delay time and the corresponding Moho depth are listed in Table 4.2. Table 4.2: Ps-P delay time and Moho depth Stations Back azimuthal interval P and Ps delay time (sec) Moho depth (km) KNN NE 4.95 40.7 SE 4.75 39.1 SW 4.95 40.7 NW 4.75 39.1 KBB NE 3.88 30.4 SE 3.85 31.7 SW 5.20 42.3 NW 4.70 38.7 92 4.3.2.2. Inversion results of average receiver function at KNN and KBB The receiver functions were inverted using software ?manyinv? based on the time domain waveform inversion scheme of Owens [1984]. We computed the final inversion models for the average receiver function obtained from NE, SE, SW and NW quadrants around KNN and KBB. Due to the non-uniqueness in the inversion of receiver functions, eight inversions were performed for each run with three different initial models. The thickness of layers in the initial model varied from 1.5 km between 0 to 1.5 km depth, 2.5 km between 1.5 to 11.5 km depth, 3.0 km between 11.5 km to 32 km depth and 4 km between 32 km to 47 km depth. We sorted out the acceptable solutions in each inversion and selected those which fit about 70% or more of the signal power. Each model produces one or more synthetic waveforms that compare well with the observed average receiver function at KNN and KBB (Fig.4.6A and B). The velocity model at each station is obtained by averaging all these solutions for the three initial models. Results are shown in Fig.4.7A and B. By assuming a Poisson?s ratio equal to 0.25, the S- velocity model can be calculated using the relation 73.1 s p V V To estimate the crustal structure in the entire Virunga area, an average velocity structure was calculated using velocity structures obtained at KNN and KBB by taking the arithmetic mean of the velocity at each layer (Figure 4.7). The gross solution is shown in Figure 4.8A. 93 In order to use the velocity model in the popular hypocenter determination program Hypocenter 3.2, we simplified this model by further combining layers with small scale high velocity with those embedded in alternating layers of high and low velocities and showed only the average velocity taking account of layers thickness. The average velocity model is shown in Figure 4.8B and listed in Table 4.3 Table 4.3: Average local velocity model in the Virunga area Vp (km/sec) Top of the layer (km) 5.30 0.0 6.50 4.0 6.96 32.0 7.67 39.0 7.92 43.0 8.03 47.0 94 Figure 4.6 Synthetic waveforms produced by the three initial models which compare well to the observed average receiver function at (A) KNN and (B) KBB 95 Figure 4.7 Acceptable velocity models obtained from the three initial models (Left) and average velocity model (Right) at (A) KNN and (B) KBB stations. Each model produces a synthetic waveform that compares well with the average receiver function. The average velocity model was obtained by averaging all the acceptable solutions for the three initial models shown at left. 96 Figure 4.8 (A) Estimate of the crustal structure for the Virunga area. It was obtained by averaging the two velocity structure obtained at KNN and KBB stations. (B) Simplified crustal structure in the Virunga area. It was obtained by combining layers with small scale high velocity with those embedded in alternating layers of high and low velocities and considering only the average velocity taking account of layers thickness. 97 4.3.3. Discussion 4.3.3.1 Lateral variation in receiver function with back azimuth At first glance at Figure 4.5A-B and Table 4.2, comparison of receiver functions (RF) along different quadrants indicates that the crust-mantle boundary beneath KNN is sharper (39-41 km) than that of KBB (30-42 km). However, the crust- mantle boundary for the SW and NW quadrant of KBB is similar to those obtained at NE and SE quadrant at KNN station. Figure 4.8C shows the Moho piercing points of rays at a depth of 36 km. The Moho piercing points were computed using IASP91 model in the taup_pierce software programme (Crotwell et al., 2001). The receiver function in the NE and the SE quadrant obtained at KBB station may also sample the Lake Kivu area and vicinity in the Western Rift. The value of 30 km for crustal thickness observed in these two quadrants is quite similar to that observed from the lithospheric cross section of the Northern Kenya rift (Achauer et al., 1992), (see also Figure 4.9). Figure 4.9 shows the gravity data obtained from the western part of the Virunga volcanic complex. It shows that the negative Bouguer gravity anomalies decrease from east to west. The two-dimensional model of these anomalies indicates that a surface layer of low-density material fills a V-shaped trough structure beneath the Lake Kivu area. KNN station is located on the western edge of the V-shaped trough. The calculated average crustal thickness beneath KNN station is 40 km. The crustal thickness values observed beneath stations KNN and KBB are in agreement with this density model shown in Figure 4.9. 98 Figure 4.8C A satellite image of the Virunga volcano group and Lake Kivu. The filed triangles and circles indicate seismographic stations and Moho piercing points of rays at 36 km depth of teleseismic events used for receiver function analysis. 99 Figure 4.9: (a) Map of Bouguer gravity anomalies (5 mgal) contour interval. Arrows= eastern limit of Precambrian basement, Triangles= location of volcanic cones. (b) Two dimensional model (dashed line) of Bouguer anomalies observed along D-D? (circles) shown above density model and topographic profile. Source: Zana et al. (1992b) 100 4.3.3.2 Crust-mantle boundary beneath KNN and KBB As stated before, these two stations are only 29 km apart, but the crustal structure beneath these stations is quite different. It is noteworthy that KNN station is close to volcano Nyamuragira. In contrast, the KBB station is near volcano Nyiragongo (See Figure 4.1 and 4.8C). However, a glance at their structures (Figure 4.7A-B.) shows that they present some similarity. A zone of high velocity anomaly is observed in both cases. It is observed at KNN from 3-20 km with an average velocity of 6.9 km/s and from 3- 10 km with an average velocity of 7.3 km/s at KBB. A low velocity zone is observed at KNN and KBB from 20-30 km and 18-28 km, respectively, with average velocities of 6.1 and 5.9 km/s. At KNN, beneath 30 km, the velocity increases gradually with depth at a constant rate of 0.46 km/s. This rate suddenly jumps to 0.60 km/s from 39-43 km. At greater depths the velocity increases again with a small rate of 0.18 km/sec and then remains constant with a value of about 7.9 km/s. At KBB, a jump in the rate of increase in velocity from 0.19 km/s to 0.46 km/s is observed from 28-30 km. At greater depths the rate decreases from 0.33 km/s to 0.05 km/sec, and then jumps again from 0.05 to 0.59 km/sec from 35-39 km. These zones between 39-43 km and 30-39 km may likely represent the crust-mantle boundary beneath KNN and KBB, respectively. From P and Ps delay times (Table 4.2), we obtained a crust-mantle transition zone with a depth of about 39-41 km and 30 km at KNN and KBB, respectively 101 4.3.3.3 High and low velocity zone A high-velocity anomaly (6.9 to 7.3 km/s) is observed beneath stations KNN and KBB at depth about 3-20 km and 3-10 km, respectively. High-velocity regions around active volcanoes are typically interpreted to indicate intrusive complexes of basaltic dikes and ultramafic rocks [Thurber, 1984; Okubo et al., 1997; Hill and Zucca; Park et al., 2007]. The average value of the high-velocity material at these stations is different because receiver functions at these stations sample quite two different structures that are influenced by the proximity to the volcano and by the material produced by each volcano. In fact, the isotopic and geochemical study of the Nyiragongo and Nyamuragira volcanics by Ramananda et al. [2007] proved that the Nyiragongo and Nyamuragira volcanics are compositionally distinct although spatially they are only 15 km apart. The Nyiragongo rocks are silica undersaturated and unusually enriched in alkalis with composition similar to foidites. In contrast, the Nyamuragira volcanics range from basalts and basanites to tephrites, similar to other volcanics around the Tanzanian Craton. To account for the major trace element as well as isotopic composition of the Nyiragongo and Nyamuragira rocks, Ramanda et al. [2007] presented a model where the Nyiragongo lavas were derived from comparatively greater depth by very low degree partial melting of a carbonate-metasomatized and phlogopite-bearing garnet peridotite source assemblage. The Nyamuragira lavas are the products of a larger degree of partial melting of the same source, but from comparatively shallow depths. This model is consistent with the different values of 6.9 km/s and 7.3 km/s of body material found around KNN and KBB, respectively. However, 102 to confirm this result, the exact information on rock velocity of volcanoes Nyiragongo and Nyamuragira is needed. These data are not available. The low velocity zone found beneath KNN and KBB at depths from 20-30 km and 18-28 km, respectively, may probably related to a magma chamber or a melt- rich sill. This interpretation is supported by the deep long-period volcanic earthquakes observed 10 or 11 months prior to the next Nyamuragira eruptions from a deep source located at about 20-30 km (Mavonga et al., 2006). 4.3.3.4 Application of the new velocity model to the hypocenter determination The new average velocity model (Figure 4.8B) was applied to the hypocenter determination of hybrid (low- and high-frequency) volcanic earthquake swarms with clear P-wave onset prior to the 27 November 2006 eruption of volcano Nyamuragira. We used the Hypocenter 3.2 location program by Lienert and Havskov [1995] and considered only hypocenters with a standard root mean square of travel time residual (rms) less than 0.3 sec, standard error in latitude and longitude less than 5 km and maximum standard error on focal depth equal to 7 km. Hypocenter maps were obtained using both the new model and the model used at Goma Volcano Observatory GVO [Mavonga et al., 2006], see Figure 4.10A-B The hypocenter maps were compared with each other, together with InSAR data provided by National Museum of Natural History, Luxembourg. The results of the emergency InSAR acquisition and processing data [D'Oreye et al., 2007] showed that there was inflation along the axis of a fracture linking the two volcanoes and a subsidence pattern on the East of Nyamuragira. This deformation extended up to 103 the south-western flank of Nyiragongo. The epicentres determined using the new velocity model show better agreement with the InSAR data and the eruption site than the epicentres provided by the GVO model. Figure 4.10 Hypocenters map of volcanic earthquakes for the period from 26 November 2006 to 27 November 2006 computed using (A) GVO velocity model and (B) the new velocity model .The open circles and lozenge mark the earthquake epicentres and eruptive site, respectively. The coordinates of the eruptive site were determined by GPS 104 In fact, the eruptive site is situated at the intersection of the fracture linking the two volcanoes and a fracture trending to the southwest flank of Nyiragongo. Furthermore, these epicentres are distributed along a linear trend, which may be associated with the opening of fissures before eruption. 4.4 Seismic activity prior to the recent eruptions of volcano Nyamuragira 4.4.1 Method In this study, only analogue records were used to locate events because KTL and LBG stations (the stations nearest to the central crater of Nyamuragira volcano) were not yet relayed to Goma base station due to transmission problems. The digital records were, however, used to estimate the predominant frequency content of the waveforms used in this study from November 2003 to November 2006 (Figure 4.11). On the basis of the waveform pattern, the volcanic events were tentatively classified as follows [Fehler and Chouet, 1982; Tanaka, 1983, McNutt, 1992 and Lukaya et al., 1992]: Short-period (SP) earthquakes: Earthquakes having P- and S-phases discernible as tectonic earthquakes and predominant high-frequency component content greater than 5 Hz. Their S-P times are less than 5 sec and are located in the Virunga volcanic area. 105 Long-period (LP) earthquakes: Transient signals having weak P- and emergent or no S-phases with a predominant low frequency content component between 1 and 3 Hz. Volcanic tremor: A sustained occurrence of long-period events appearing on a seismogram as an irregular sinusoid of long duration compared with a tectonic earthquake of the same amplitude. We used the Hypocenter 3.2 location program by Lienert and Havskov [1995] to locate events. The periods of observation were 18 August 2002 to 7 May 2004 and 1 July 2004 to 27 November 2006. The accuracy of the phase reading in the analogue records was about ?0.1 sec and ?0.5 sec for short and long-period events, respectively. In this study, we considered a combination of the crust model obtained from a simple two-layer model that is the average of the Bonjer et al., [1970], Bram [1975] and Nolet and Mueller [1982] models for the Western Rift Valley of Africa, and that obtained from teleseismic P-wave receiver function analysis beneath two seismic broadband stations in the Virunga area [this study]. The P- wave velocity is 4 km/s in the uppermost 3 km of our model. From 3-4 km and 4- 20 km, we adopt a velocity of 4.49 km/sec and 6.35 km/sec, respectively. From 20 km to the Moho (~ 30 km), we specify a P-wave velocity of 6.65 km/sec. For the upper mantle (> ~30 km), a P-wave velocity of 6.9 km/sec and 6.96 km/sec are given from 30-32 km and 32-39 km, respectively. The S-wave velocities are calculated using the Vp/Vs ratio of 1.73 (Figure 4.12). 106 (a) (b) Figure 4.11 Velocity spectra of typical (a) LP events (18 March 2004) and (b) hybrid (LP + short) events (27 November 2006) 107 Figure 4.12: Velocity models used in this study. ABBN-MOD.vp = Average model of Bonjer et al. (1970), Bram (1975) and Nolet and Mueller (1982) for the Western Rift Valley of Africa; RF-MOD.vp= velocity model obtained in this study from teleseismic P-wave receiver function analysis; combinmod.vp= combination model of both these models. 108 To reduce bias due to uncertainty in the phase reading and velocity model, we considered only epicentres with standard error in latitude and longitude less than 3.8 km and a standard root mean square (rms) error on the travel-time residual less than 0.15 sec. As our source region is located inside the network, the epicentre determinations were generally good. The focal depths are less accurate, especially for deeper earthquakes (in the range of 10-30 km). We adopted an optimal maximum standard error in the focal depth erz=7 km (see Mavonga et al. 2006). 4.4.2. Temporal and spatial variations in seismicity Figure 4.13A-B shows the monthly count of locatable (i.e. registered by 3 or more stations) and located volcanic earthquakes in the eleven months prior to the 8 May 2004 and 27 November 2006 eruptions of Nyamuragira. These figures show that these eruptions were preceded by a progressive increase in the number of long- period earthquakes over an eleven month period. Figures 4.14 to 4.16 show the three stages of seismicity activity leading up to 2004 and 2006 eruptions. 1. The activity at stage 1 (Figure 4.14) is probably related to the post-eruptive activity associated with the previous eruption of Nyamuragira. Seismicity was dominated by the occurrence of shallow (0-5 km) long-period (LP) events. 2. During the second stage (Figure 4.15), which commenced 10 or 11 months prior to the next eruption, swarms composed mainly of LP events were observed. Some local tectonic earthquakes were also felt by inhabitants in the Rutsuru and Virunga areas. Seismicity was dominated by the presence of shallow (0-5 km) and deep (20-30 km) LP events. An aseismic zone was observed between the shallow and deep cluster at depths of 4-8 km. This aseismic zone is probably a region of low rigidity occupied partly by magma, 109 similar to that detected by Hamaguchi et al. [1983]. The deep LP volcanic earthquakes observed at this stage can be attributed to a fairly steady flow of magma into the Nyamuragira shallow reservoir from a deep source located at about 20-30 km. This deep zone coincides with the low-velocity zone observed beneath KNN station [This study]. 3. In the final stage (Figure 4.16), 1 or 2 months before the next eruption, most LP events clustered at a shallow depth. These events are probably associated with the filling of the shallow reservoir with magma to its threshold volume. 4. 4.3 Summary of Observations Eleven months before the 8 May 2004 and 27 November 2006 eruptions of Nyamuragira, an increase in seismicity related to the number of long period was observed. This increase was due to swarm-type seismicity composed mainly of long?period events observed during this period. Shallow long-period earthquakes, typifying repose-period seismicity (2002-2003 and 2004-2005), were gradually augmented by deep long- period events ten or eleven months prior to the eruption. The number of deep long-period events declined in the final month prior to the eruption, while the number of shallow long-period events increased. These two aspects of seismicity of long-period events (i.e. changes in the frequency of occurrence and in the depth), integrated with other available data (e.g. INSAR, GPS, geochemical, geological) can be used to characterize volcanic processes and forecast volcanic eruption in the Virunga area. 110 Figure.4.13: Monthly count of locatable (series 1) and located (series 2) long- period volcanic earthquakes 11 months before the 8 May 2004 eruption of Nyamuragira (Top) and the 27 November 2006 eruption of Nyamuragira (Bottom) . 111 Figure 4.14: Long and short-period earthquake hypocenters observed in the Virunga area during the period (A) 24 November 2002 to 24 December 2002 and (B) 01 July 2004 to 01 August 2004. The circle defines the area used for the depth section within which epicentres are selected. Epicentres within the circle are projected normal onto the cross-section. The continuous blue line crossing the epicentres is the line of cross-section of earthquake hypocenters. It is defined by a centre marked by a large star (set at the summit of the Nyamuragira crater) and an azimuth (set to 160?).This direction coincides with that of the main fissure connecting the volcanoes Nyiragongo and Nyamuragira. The large filled triangles indicate seismic stations. The small filled red circles and green squares mark long and short- period events, respectively. eismicity was dominated by shallow (0-5km) LP events. This activity is observed immediately after the previous eruption. 112 Figure 4.15: Long and short-period earthquake hypocentres observed in the Virunga area during the period (A) 07 January 2004 to 07 February 2004 and (B) 25 July 2006 to 25 August 2006. Seismicity was dominated by shallow and deep LP events (0-5 km and 20-30 km). This activity is observed about 10 months before the next eruption. 113 Figure 4.16: Long and short-period earthquake hypocentres observed in the Virunga area during the period (A) 07 April 2004 to 07 May 2004 and (B) 28 October 2006 to 27 November 2006. A large increase in shallow LP events and a reduction in deep LP events were observed. This activity is observed about 1 month before the next eruption 114 CHAPTER 5 CONCLUSIONS AND RECOMMENDATIONS 5.1. Conclusions From the probabilistic seismic hazard analysis in DR Congo and surrounding areas, four seismic hazard source zones have been identified and rated as follows: 1. Zone A (very high hazard), largely the Lake Tanganyika Rift zone where the PGA values of 0.32g, 0.22g and 0.16g are expected to occur with probability 2%, 5%, and 10 % in 50 years, respectively. 2. Zone B (high hazard), which includes the Lake Kivu basin, Ruwenzori and Lake Edouard region. 3. Zone C (moderate hazard), which includes Rutsuru, Masisi, Upemba and a part of the Congo basin close to Western Rift . 4. Zone D (low hazard), which includes the remainder of the Congo basin.. This division is based on assessment of the seismicity and the expected intensity of ground motion. Each seismic hazard zone can be defined by a seismic hazard zoning coefficient. Usually the strongest earthquake zone coefficient is assumed to be 1.0. Two studies were conducted to understand how volcanoes work and reduce the risk posed by the Virunga volcanic eruptions in the Western Rift Valley of Africa: 1. Firstly, we estimated the crustal structure beneath KNN and KBB broadband stations located in the Virunga volcanic area by inverting stacks of teleseismic receiver functions. From this study, it was revealed: 115 A lateral variation in receiver function with back azimuth, implying lateral heterogeneities or dipping of layers, was observed. The crust-mantle transition zone beneath the Virunga area is estimated at a depth from 30 to 43 km. A low-velocity zone was observed around KNN and KBB at depth from 20 to 30 km and 18 to 28 km, respectively. This zone may probably relate to a magma chamber or a melt-rich sill in the field of Nyamuragira and Nyiragongo volcanoes. A high-velocity zone with an average velocity of 6.9 km/sec was found at a depth of 3 to 20 km beneath KNN station. A similar high-velocity zone with an average velocity of 7.3 km/sec at a depth of 3 to 10 km was also found at KBB. These zones are indicative of magma cumulate within the volcanic edifice. 2. Using the new average velocity model in the Virunga area, the spatial and temporal variations in the seismicity in the Nyamuragira area were examined for the periods 18 August 2002 to 7 May 2004 and 1 July 2004 to 27 November 2006, prior to the 8 May 2004 and 27 November 2006 eruptions of Nyamuragira, respectively. It was found that the seismicity exhibits similar tendencies during both periods: Swarm-type seismicity, composed mainly of long period earthquakes, preceded both eruptions of Nyamuragira and was probably enhanced by tectonic seismicity related to rifting. 116 Ten or eleven months before eruption, a steady increase in seismicity at a constant rate from a deep magma feeder was observed. In the last stage (1 or 2 months) before the eruption, the hypocenters of long-period earthquakes became shallower. 5.2. Recommendations The source zonation can be improved by supplementing the area sources used in this study with fault sources. An effective earthquake disaster mitigation strategy requires that base maps of known faults be compiled, and efforts to detect possible unknown faults be made using paleoseismology. Paleoseismology is the science to study the nature, timing and location of pre-instrumental earthquakes. If paleoseismic events are well documented, we can evaluate earthquake potential of specific faults. We need information on paleo-earthquakes for a sufficiently long period of time because active faults on land have long recurrence cycles. Evidence of paleo-earthquakes can be found by investigating geologic sequences (depositional sediments) through direct excavation survey across active faults. The PGA exceeding indicated values determined in this study are only controlled by magnitude, distance and regional variation of attenuation of strong motion. The local geology (soil type) was not taken account. As we know, the near-surface geological condition can control the amplification of ground motion at a particular site (Jongmans, 1989; Jennings, 1971; 117 Poccski, 1969; Sozen, 1968; Celebi, 1987). Damage patterns in past earthquakes show that soil conditions at a site may have a major effect on the level of ground shaking. The characteristics of ground surface motion (i.e PGA) are strongly dependent on the local soil conditions. Therefore, a local site effect analysis should be conducted to evaluate the amplification of seismic wave ground motion by soft near-surface deposit to the underlying bedrock, or firm soil considered as rock. Mapping of seismic hazard at local scales to incorporate the effects of local soil conditions is called microzonation for seismic hazard. The next step will be the development of a frequency-dependent attenuation relationship. The result of seismic hazard analysis should be convolved with a seismic fragility function, which quantifies the probability of various level of damage to a facility as a function of ground motion. This will give a seismic risk analysis which indicate probabilities per unit time of different levels of failure or loss The average velocity structure beneath Virunga volcanic area obtained in this study can be improved by setting up two supplementary broadband stations at RSY and KTL and using a long study period. It can be used as a first approximation for high resolution tomography in the Virunga area. 118 REFERENCES Abramowitz, M., Stegun, J.A., 1970. Handbook of Mathematical Functions (9th ed.), Dover Publ., New York. Achauer, U., Maguire, P., Mechie, J., Green, W.V., The KRISP Working group, 1992. Some remarks on the structure and geodynamics of the Kenya Rift. Tectonophysics 213, 257-268. Ambraseys, N.N., Adams, R.D., 1992. Reappraisal of major African earthquakes, south of 20?N, 1900-1930. Tectonophysics, 209, 293-296 Ammon, C.J., Randall, G.E., Zandt, G., 1990. On the nonuniqueness of receiver function inversions. J. Geophys. Res. 95, 15303-15318. Ammon, C.J., 1991. The isolation of receiver effects from teleseismic P-wave fronts. Bull. Seismol. Soc. Am. 81, 2504-2510. Ammon, C.J., Zandt, G., 1993. Receiver structure beneath the southern Mojave block, California. Bull. Seismol. Soc. Am., 83, 737-753. Ammon, C.J., 1997. Research notes and softwares. Available at http://eqseis.psu.edu/ cammon/ Anderson, J. A., Wood, H. O., 1925. Description and theory of the torsion seismometer, Bull. Seism. Soc. Am. 15, 1-72. Asfaw, L., Bilham, R., Jackson, M., Mohr, P., 1992. Recent inactivity in the African Rift, Nature, 357, 447. Asish, R.B., Aba, P.S., Ramananda, C., Orlando, V., Dario, T., 2007. Petro- chemical geodynamics of Nyiragongo and Nyamuragira volcanoes in the Western Rift of the East African Rift System. Proceeding of the 26th ECGS workshop, AVCOR07, Luxemburg, November 2007. 119 Atalay, A., 2002. Active compressional tectonics in Central Africa and implications for plate tectonic models: Evidence from fault mechanism studies of the 1998 earthquakes in the Congo Basin, Journal of African Earth Sciences, 35, 45-50. Atkinson, G.M., Boore, D.M., 2006. Earthquake ground motion prediction equations for Eastern North America, Bull. Seism. Soc. Am., 96, 2181- 2205. B?th, M., 1975. Seismicity of the Tanzania region. Tectonophysics, 27, 353-379. B?th, M., 1981. Earthquake magnitude recent research and current trends. Earth Sci. Rev., 17, 171- 186. Bender, B., Perkins, D. M., 1987. Seisrisk III-A computer program for seismic hazard estimation, US Geological Survey Bulletin, 1772, 48 pp. Bilham, R., Bendick, R., Larson, K., Mohr, P., Braun, J., Tesfaye, S., Asfaw, L., 1999. Secular and tidal strain across the Main Ethiopian Rift, Geophysical Research Letters 26, 2789-2792. Bonjer,K.P.,Fuchs,K., Wohlenberg,J., 1970 .Crustal structure of the East African Rift System for spectral response Ratios of long-period body waves, Zeitshr Geophys., 36, 287-297. Boore, D.M., 1983. Stochastic simulation of high frequency ground motions based on seismological models of the radiated spectra, Bull.Seism.Soc.Am. 73, 1865-1894. Boore, M. D., Joyner, W. B., 1984. A note on the use of random vibration theory to predict peak amplitudes of transcient signals, Bull. Seism. Soc. Am. 74, 2035-2039. 120 Bram, K., 1972. Seismicity of Katanga and Western Zambia, Southwest of East African Rift System, from 1960 to 1971, Bull. Seism. Soc. Am. 62, 1211- 1216. Bram, K., 1975. Zum Aufbau der Kruste und oberen mantels in bereich des weslchen grabens des Ostrafrikansichen Grabensystem und im Ostlinchen Zaire-Bechen. Ergebnisse Einer Untersuchung der aumwellem von Nah-Erdbeben, Geophysikalische Abhandlungen , Heft.4, 1-65. Cahen, L., 1952. Esquisse Tectonique du Congo Belge et Ruanda-Urundi, 1:3,000,000. Ministere des Colonies, Commission de Geologie, Bruxelles Cahen, L., 1954. Geologie du Congo Belge, Liege. Camelbeek, T., Iranga, M. D., 1996. Deep crustal earthquakes and active faults along the Rukwa trough, eastern Africa, Geophys. J. Int., 124, 612-630. Cassidy, J.F., 1992. Numerical experiments in broadband receiver function analysis. Bull. Seismol. Soc. Am. 82, 1453-1474. Celebi, M., Prince, J., Dietel, C., Onate, M., Chavez, G., 1987. The culprit in . Mexico City- Amplification of motions, Earthquake Spectra, 3, 315-328. Chiarabba, C., Amato, A., Boschi, E., Barberi, F., 2000. Recent seismicity and tomographic modeling of the Mount Etna plumbing system, J. Geophys. Res. 105: 10923-38. Chorowicz, J., 2005. The East African rift system, Journal of African Earth Sciences 43, 379-410. 121 Chu, D., Gordon, R.G., 1999. Evidence for motion between Nubia and Somalia along the Southwest Indian ridge. Nature 398, 64-66. Clayton, R.W., Wiggins, R.A., 1976. Source shape estimation and deconvolution of teleseismic body waves. Geophys. J. R. Astrom. Soc. 47, 151-177. Clay, C., Cassandra, R.C., 1995. Volcanism in Eastern Africa, Final report, NASA/ASEE Summer Faculty Fellowship Program, Johnson Space Center, Contract Number: NGT-44-001-800 Cornell, C. A., 1968. Engineering seismic risk analysis, Bull. Seism. Soc. Am., 18, 1583-1606. Cramer, H., 1961. Mathematical methods of statistics. Princeton University Press, Princeton. Crotwell, H.P., Owens, T.J., Ritsema, J., 2001. Flexible seismic travel time and ray path utilities, User?s guide. Department of Geological Sciences, University of Carolina. Available at http://www.seis.sc.edu Daly, M. C., Chorowicz, J., Fairhead, J. D., 1989. Rift basin evolution in Africa: the influence of reactivated steep basement shear zones, Cooper, M. A. and Williams, G. D. (Eds.): Inversion tectonics, Geological Society, London, Special Publication, 44, 309-335. Daly, M.C., Lawrence, R., Kimuna, D., Benga, M., 1992. Tectonic evolution of the Cuvette centrale, Zaire. J. Geol. Soc., London, 149, 539-546 De Bremaecker, J. Cl., 1955. Determination des magnitudes des seismes du Congo Belge, Acad. Roy. Soc. Col. Bull. Seances, 1, 1043-1046. De Bremaecker, J. Cl., 1959. Seismicity of the Western African Rift Valley, J. Geophys. Res. 64, 1961-1966. 122 Delvaux, D., 2001. Karoo rifting in Western Tanzania precursor of Gondwana break-up, Contribution to the Geology and Paleontology of Gondwana, in: Honour of Helmut Wopfner, Cologne, 111-125. Demant, A., Lestrade, P., Lubala, R.T., Kampunzu, A.B., Durieux, J., 1994. Volcanological and petrological evolution of Nyiragongo volcano. Virunga volcanic field, Zaire. Bull. Volcanol., 56, 47-61 De Mets, C., Gordon, R.G., Argus, D.F., Stein, S., 1990. Current plate motions, Geophysical Journal International, 101, 425-478. D?Oreye, N., Cayol, V., Kervyn, F., GVO team, 2007. The November 2006 Nyamulagira eruption revealed by lnSAR . Proceeding of the 26th ECGS workshop, AVCOR07, Luxemburg, November 2007. Douglas, J., 2007. On the regional dependence of earthquake response spectra, ISET Journal of Earthquake Technology, Paper No 477, Vol. 44, No 1, pp. 71-99. Dugda, M.T., Nyblade, A.A., Julia, J., Langston, C.A., Ammon, C.J., Simiyu, S., 2005. Crustal structure in Ethiopia and Kenya from receiver function analysis: Implications for rift development in eastern Africa. Journal of Geophysical Research, 110, B01303, doi:10.1029/2004JB003065. Dunbar, J. A., Sawyer, D. S., 1989. How pre-existing weakness control the style of continental breakup, Journal of Geophysical Research, 94, B6, 7278- 7292. Dziewonski, A. M., Elstrom, G. Salgamik, M. P., 1996. Centroid moment tensor solutions for July- September 1995, Phys. Earth Planet. Interiors 97, 3- 13. Ebinger, C. J., 1989a. Tectonic development of the Western branch of the East 123 African Rift system, Geol. Soc. of Am. Bull., 101: 885-993. Ebinger, C. J., 1989b. Geometric and kinematics development of border faults and accommodation zones, Kivu - Ruzizi rift, Tectonics, 8: 117-133. Fairhead, J. D., Girdler, R. W., 1972. The seismicity of the East African rift system, Tectonophysics, 41,19-26. Fairhead, J. D., Henderson, N. B., 1977. The seismicity ofsouthern Africa and incipient rifting, Tectonophysics, 41, 19-26. Fairhead, J. D.,Stuart, G. W., 1985. Seismicity of the East Africa Rift system and comparison with other continental rifts. In:Palmason, G. (Ed.), Continental and oceanic rifts, geodynamic series, vol.8, pp. 41-61. Fehler, M., Chouet, B., 1982. Operation of a digital seismic network on Mount St Helens volcano and observations of long-period seismic events that originate under the volcano, Geophys. Res. Letter, 9, 1017-1020. Ferdinand, R.W., 2007. Seismic hazard analysis of Arusha , Tanzania. Abstract, AfricaArray Workshop, 17-18 July 2007, Johannesburg, South Africa Foster, A.N., Jackson, J.A., 1998. Source parameters of large African earthquakes: implications for crustal rhelogy and regional kinematics. Geophysical Journal International, 134, 422-448. Frankel, A., 1995.Simulating strong motions of large earthquakes using recordings of small earthquakes, Bull. Seism. Soc. Am. 85, 1144-1160. Frohlich, C., Davis, S. D., 1993. Teleseismic b values or, much ado about 1.0, J. Geophys. Res., 98, 631-644. Furon, R., Daumain, G., 1959. Esquisse structural proviso ire de l?Afrique: 1:10-106 , Congr. Geol. Intern. Assoc. Serv. Geol. Africains, Paris. 124 Girdler, R.W., McConnell, D.A., 1994. The 1990 to 1991 Sudan earthquake sequence and the extent of the East African rift system, Science, 264 : 67-70. Gumper, F. and Pomeroy, P.W., 1970. Seismic wave velocities and earth structure in the African continent. Bull. Seism. Soc. Am. 60,651-668. Gutenberg, B., Richter, C.F., 1949. Seismicity of the Earth, Princeton University Press, Princeton, New Jersey, USA. Gutenberg, B., Richter, C. F., 1954. Seismicity of the Earth and associated phenomena, Princeton University Press, New Jersey. Gurrola, H., Minster, J.B., Owens, T., 1994. The use of velocity spectrum for stacking receiver functions and imaging upper mantle discontinuities. Geophys. J. Int. 117, 427-440. Hamaguchi, H., 1983. Seismological evidence for magma intrusion during 1981- 1982 Nyamuragira eruption. In: Hamaguchi, H.( Editor), Volcanoes Nyiragongo and Nyamuragira geophysical aspects, Tohoku University, Sendai, Japan, pp35-42. Hanks, T., Kanamori, H., 1979. A moment magnitude scale, J. Geophys Res., 84, 2348-2350. Hartnady, C.J.H., Benouar, D., 2007. African catalogue of earthquakes (ACE) project : towards earthquake risk reduction in active plate-boundary zones, AfricaArray Workshop, 17-18 July 2007, Abstracts. Hartzell, S.H., 1978.Earthquake aftershocks as Green's functions, Geophys. Res. Lett. 5, 1-4. Havskov, J., Ottemeller, L., 2006. Seisan: The earthquake Analysis software 125 (Version 8.1), User manual. Hayashi, S., Kasahara, M., Tanaka, K., Hamaguchi, H., 1992 .Major elements chemistry of recent eruptive products from Nyamuragira volcano, Africa (1976-1989). In: Hamaguchi, H. (Editor), Geophysical study on the hotspot volcanoes in the African Continent, Publ. Faculty of .Science, Tohoku University, Sendai, Japan, 83-87. Hill, D.P., Zucca, J.J., 1987. Geophysical constraints on the structure of the Kilauea and Mauna Loa volcanoes and some implications for seismomagmatic processes, U.S. Geol. Surv. Prof. Pap. 1350, 903-917. Hlatywayo, D. J., 1992. Seismicity of Zimbabwe and surrounding areas during the period 1959-1990, Seismol. Dept., Uppsala, Rep., 3-92, 1-98. Hlatywayo, D. J., 1995. Seismotectonics and seismic hazard estimates in Central and Southern Africa. Acta Univ. Ups., Comprehensive summaries of Uppsala dissertation from the Faculty of Science and Technology 160, 12pp, Uppsala. Irikura, K., 1983. Semi-empirical estimation of strong ground motions during largeearthquakes, Bull. Disaster Prevention Res. Inst., Kyoto Univ., 33, 63-104. Irikura, K., 1986. Prediction of strong acceleration motions using empirical Green'sfunction, Proc.7th Japan Conf. Earthquake Engineering, 151- 156. Jennings, P.C., 1971. San Fernando earthquake of February 9, 1971, EERL, 71- 02, Calif. Inst. of Techn., Pasadena. Jestin, F., Huchon, P., Gaulier, J.M., 1994. The Somalia Plate and East African 126 Rift system: present-day kinematics. Geophysical Journal International 116, 637-654. Jonathan, E., 1996. Some aspects of seismicity in Zimbabwe and eastern and southern Africa. Msc. Thesis, Institute of Solid Earth Physics, Univ.Bergen, Norway, 100 pp. Jongmans, D., 1989. Les phenomenes d?amplification d?ondes sismiques dus a des structures geologiques, Annales de la Societe Geologique de Belgique, Tome 112 (fascicule 2), pp 369-379. Joyner, B. W., 1984. A scaling law for the spectra of large earthquakes. Bull. Seism. Soc. Am., 74, 1167-1188. Kampunzu, R. B., Caron, J. P. H., Lubala, T. R., 1986. The East African rift, magma genesis and astheno-lithospheric dynamics. Episodes 9, 211-216. Kaspar, T., Ritter, J.R.R., 1998. P-SV conversions of teleseismic waves beneath Chyulu hills volcanic field, Kenya, Geophys. Res. Let. 25, 559-562. Kavotha, K.S., Mavonga, T., Jacques, D., Mukambilwa, K., 2003. Towards a more seismic picture of the January 17th, 2002 Nyiragonga eruption. Acta Vulcanologica Vol. 15 (1-2), 87-100. Kijko, A., Sellevoll, M. A., 1989. Estimation of earthquake hazard parameters from incomplete data files. Part I: Utilization of extreme and complete catalogs with different threshold magnitudes, Bull. Seism. Soc.Am.,79, 3, 645-654. Kijko, A., Graham, G., 1998. Parametric-historic procedure for probabilistic seismic hazard analysis. Part I: Assessment of maximum regional magnitude Mmax, Pure appl. Geophys. 152, 413-442. 127 Kijko, A., Graham, G., 1999. Parametric-historic procedure for probabilistic seismic hazard analysis, Part II : Assessment of seismic hazard at specified site, Pure Appl. Geophys., 154, 1-22. Kijko, A., Retief, S.J.P., Graham, G., 2002. Seismic hazard and risk assessment for Tulbang, South Africa: Part I- Assessment of seismic hazard, Natural Hazards, 26, 175-201. Kijko, A., Graham, G., Bejaichund, M., Robin, D., 2003. A probabilistic seismic hazard map for Sub-Saharan Africa, Report Number: 2003-0035, Council for Geoscience, South Africa. Kijko, A., 2004. Estimation of the Maximum earthquake magnitude Mmax, Pure appl. Geophys. 161, 1655-1681. Kinabo, B.D., Atekwana, E.A., Hogan, J.P., Modisi, M.P., Wheaton, D.D., Kampunzu, A.B., 2007. Early structural development of the Okavango rift zone, NW Botswana, Journal of African Earth Sciences, Volume 48, Issues 2-3, 125-136. Kind, R., Vinnik, L.P., 1988. The upper mantle discontinuities underneath the GRF array from P to S converted phases, J. Geophys. 62, 138-147. Knopoff, L., Kagan, Y. Y., 1977. Analysis of the theory of extreme as applied to earthquake problem. J. Geophys. Res. 82: 5647-5657. Komorowski ,J.C.,Tedesco, D., Kasereka, M., Allard, P., Papale, P., Vaselli, O., Durieux, J., Baxter, P., Halbwachs, M., Akumbi, M., Baluku, B., Briole, P., Ciraba, M., Dupin, J.C., Etoy, O., Garcin, D., Hamaguchi, H., Houlie, N., Kavotha, K.S., Lemarchand, A., Lockwood, J., Lukaya, N., Mavonga, T.G., De Michele, M., Mpore, S., Mukambilwa, K., Munyololo, F., 128 Newhall, C., Ruch, J., Yalire, M., Wafula, M., 2003. The January 2002 flank eruption of Nyiragongo volcano (Democratic Republic of Congo) : chronology, evidence for a tectonic rift trigger, and impact of lava flows on the city of Goma, Acta Vulcanologica, Vol. 15 (1-2), 27-42. Krinitzky, E.L., Chang, F.K. and Nuttli, O.W., 1988.Magnitude related earth- earthquake ground motion. Bull. Ass. Eng. Geol. Vol. XXV, 399-423. Kubanza, M., Nishimura, T., Sato, H., 2006. Spatial variation of lithospheric heterogeneity on the Globe as revealed from transverse amplitudes of short-period teleseismic P-waves, Earth Planets Space, 58, e45-e48. Langston, C.A., 1979. Structure under Mount Rainier, Washington, inferred from teleseismic body waves, J. Geophys. Res. 84, 4749-4762. Langston, C.A., 1989. Scattering of teleseismic body waves under Pasadena, California. J. Geophys. Res. 84, 4749-4762. Lienert, B.R.E, Havskov, J., 1995. A computer program for locating earthquakes both locally and globally, Seismological Research Letters, 66, 26-36. Lombe, K. D., Mubu, M. S., 1992. Instrumentation and seismicity in Zambia, Tectonophysics, 209, 31-33. Lowrie, W., 1997. Fundamentals of Geophysics. Cambridge, pp. 14-15. Lukaya, N., Ciraba, M., Mavonga, T., Wafula, M., 1992. Main pattern of waveforms observed in the Virunga volcanic zone, Western Rift Valley of Africa, Tectonophysics, 209, 261-265. Maasha, N., Molnar, P., 1972. Earthquake fault parameters and tectonics in Africa. J. Geophys. Res.,77, 5731-5743. Maasha, N., 1975. The seismicity of the Ruwenzori region in Uganda, J. 129 Geophys. Res., 80, 1485-1496. Mangino, S.G., Zandt, G., Ammon, C.J., 1993. The receiver structure beneath Mina, Nevada. Bull. Seismol. Soc. Am. 83, 542-560. Marshall, P. D., 1970. Aspects of spectral differences between earthquakes and underground explosions, Geophys. J. R. Astr. Soc. 20, 397-416. Mavonga, T., Sadaka, K.K., Nyombo, L., Osodundu, E., Jacques, D., 2006. Seismic activity prior to the May 8, 2004 eruption of volcano Nyamu- ragira, Western Rift Valley of Africa. Journal of Volcanology and Geothermal Research 158, 355-360. Mavonga T., 2007a. Some characteristics of aftershock sequences of major earthquakes from 1994 to 2002 in the Kivu Province, Western Rift Valley of Africa. Tectonophysics, vol. 439, pp. 1-12. Mavonga, T., 2007b. An estimate of the attenuation relationship for the strong ground motion in the Kivu Province, Western Rift Valley of Africa, Phys. Earth Planet. Interiors, vol.162, pp.13-21. Mc Guire, R. K., 1976. Evaluation of earthquake risk to site, US Geol. Surv., open-file Rep., 76-77. Mc Guire, R. K., 1993. Computation of seismic hazard, Ann. Di Geofisica, 36, 181-200. McNutt, S. R., 1992. Volcanic tremor, Encyclopedia of Earth System Sciences, Vol.4, Academic Press. Midzi, V., Hlatywayo, D., Chapola, L. S., Kebede, F., Atakens, K., Lombe, D. K.,Turyomugyendo, G., Tugume, F. A., 1999. Seismic hazard assessment in Eastern and Southern Africa, Annali di Geofisica, 42, 6, 1067-1083. 130 Midzi, V., Ottemoller, L., 2001. Receiver function structure beneath three southern African broadband stations, Tectonophysics, 339, 443-454. Morley, C. K., 1988. Variable extension in Lake Tanganyika, Tectonics, 7, 4, 785-801. Mondeguer, A., Ravenne, C., Masse, P., Tiercelin, J.J., 1989. Sedimentary basins in an extension and strike slip background. The ?South Tanganyika troughs complex?, East African Rift, Bull. Soc. Geol. France, 8, V, 501- 522. Mortelmans, G., 1951. Considerartions sur la stratigraphie des terrains pre-Karroo du Nord-Ouest du Katanga , Compt. Rend. 3e Congr. Nat. Sc. 1950, 8, 26-29. Mortelmans, G., 1953. Les antecedents tectoniques du Graben de l?Upemba ( Katanga), Bull. Volcanol. Ser. 2 13, 93-98. Munyololo, Y., Wafula, M., Kasereka, M., Ciraba, M., Mukambilwa, K., Mavonga, T., Cirimwami, M., Muhagirwa, B., Bagalwa, R., Mundala, M., 1999. Recrudescence des glissements de terrain suite ? la r?activation sismique du bassin du lac Kivu,R?gion de Bukavu (Rep. Dem. Congo), Mus. roy. Afr. Centr., Dept. Geol. Min, Rapp. Ann. 1997&1998, 285-298. National Earthquake Disaster Committee, 1994. Preliminary report on earthquake disaster in Kabalore, Bundibugyo and Kasese districts. Ministry of labour and affairs, Republic of Uganda, 56p. Nolet, G., Mueller, St., 1982. A model for the deep structure of the East African Rift System from simultaneous inversion of teleseismic data, Tectono- 131 physics, 84, 151-178. Nordquist, J. M., 1945. Theory of the largest values applied to earthquakes magnitude, Trans. Amer. Geophys. Union, 26, 29-31. Okubo, P.G., Benz, H.M., Chouet, B.A., 1997. Imaging the crustal magma sources beneath Mauna Loa and Kilauea volcanoes, Hawaii. Geology 25, 867-870. Owens, T.J., 1984. Determination of crustal and upper mantle structure from analysis of broadband teleseismic P-waveforms. PhD thesis, Depatment of Geology and Geophysics, The University of Utah. Owens, T.J., Taylor, S.R., Zandt, G., 1987. Crustal structure of regional seismic test network station determined from inversion of broadband teleseismic P waveforms. Bull. Seismol. Soc. Am., 77, 631-662. Page, R., 1968. Aftershocks and microaftershocks, Bull. Seismol. Soc. Am., 58, 1131-1168. Park, J., Morgan, J.K., Zelt, C.A., Okubo, P.G., Peters, L., Benesh, N., 2007. Comparative velocity structure of active Hawaiian volcanoes from 3-D onshore-offshore seismic tomography. Earth and Planetary Sciences Letters, 259, 500-516. Pocsski, A.C., 1969. The ground effect of the Skopje July 26, 1963 earthquake, Bull. Seism. Soc. Am. 59, 1-29. Prodehl, C., Keller, G.R., Khan, M.A. (Eds.), 1994. Crustal and upper mantle structure of the Kenya Rift. Tectonophysics, 236, Special issue, 483 pp. Prodehl, C., Mechie, J., 1991. Crustal thinning in relationship to the evolution of the Afro-Arabian rift system- a review of seismic-refraction data. In: 132 Makris, J.P., Mohr and R. Rihm, (Editors), Red Sea: Birth and early history of a new oceanic basin. Tectonophysics, 198: 311- 327. Ramananda, C., Asish, R.B., Alba, P.S., Dario, T., Kenneth, W.W.S., 2007. Isotopic and geochemical study of the Nyiragongo and Nyamuragira volcanics in the Western rift, East African rift system. Proceeding of the 26th ECGS workshop, AVCOR07, November 2007. Reuter, L., 1990. Earthquake Hazard analysis: Issues and Insights, Columbia Univ. Press, New York. Richard, W.F., 2007. Seismic hazard analysis of Arusha, Tanzania, AfricaArray Workshop 17-18 July 2007, Abstracts. Richter, C. F., 1935. An instrumental earthquake magnitude scale, Bull. Seism. Soc. Am., 25, 1-32. Richter, C. F., 1958. Elementary seismology, Freeman and Co., San Francisco, 364-367. Risk Engineering, 2007. Ez-Frisk user manual (version 7.24). Risk Engineering, Inc., Golden, Colorado, USA. Risk Engineering, 2008. FRISK88M SOFTWARE user manual(version 2.18.0), Risk engineering, Inc., Golden, Colorado, USA. Rosendahl, B.R., Kilembe, E., Kaczmarick, K., 1992. Comparison of the Tanga- nyika, Malawi, Rukwa and Turkana rift zones from analyses of seismic reflection data. Tectonophysics, 213, 235-256. Sadaka, K.K., Tuluka, M., Jacques, D., Kibuye, M., 2003. Toward a more detailled seismic picture of the January 17th, 2002 Nyiragongo Eruption, Acta VulcanoLogica , Vol. 15, pp. 87- 100. 133 Sandvol, E., Seber, D., Calvert, A., Barazangi, M., 1998. Grid search modeling of receiver functions: Implications for crustal structure in the Middle East and North Africa. J. Geophys. Res. 103, 26899-26917. Sebagenzi, M. N., Vasseur, G., Louis, P., 1993. First heat flow density determina- tion from south eastern Zaire (Central Africa), Journal of African Earth Sciences, 16, 4, 413-423. Sebagenzi, M. N., Vasseur, G., Louis, P., 1997. Recent warming in South eastern Zaire (Central Africa) inferred from disturbed geothermal gradients, Paleoclimatology,Palaeoecology, Paleogeography (Global Planetary Section), 98, 209-217. Sebagenzi, M. N., Kaputo, K., 2002. Geophysical evidences of continental break up in the southeast of the Democratic Republic of Congo and Zambia (Central Africa), EGU Stephan Mueller Special Publication Series, 2, 193-206, 2002. Shudofsky, G.N., 1985. Source mechanisms and focal depths of East African earthquakes using Rayleigh-wave inversion and body-wave modelling, Geophysical Journal Astronomical Society, 83, 563-614. Singh, K. S., Ordaz, M., Lindholm, C. D., Havskov, J., 1990. Seismic hazard in Southern Norway, Bergen Univ. Seismo Series 46, 1-33. Sozen, M. A., Jennings, P.C., Mattiesen, R.B., Housner, G.W., Newmark, N. M., 1968. Engineering report of the Caracas earthquake of July 29, 1967, National Academy of Sciences, Washington DC. Somerville, P., Collins, N., Abrahamson,N., Graves, R., Saikia, C.,2001. Ground motion attenuation relations for the central and eastern United States, USGS Report No 99HQR0098, June 30, 2001. 134 Stephen, R.M., 2005. Volcanic seismology, Annual Review of Earth and Planetary Sciences, 33, 461-492. Studt, F. E., Cornet, J., Buttgenbach, H., 1908. Carte geologique du Katanga et notes descriptive, Am. Musee Congo, Ser. 2,1. Sutton, G. H., Berg, E., 1958. Seismological studies of the Western Rift Valley of Africa, Trans. Am. Geophys. Union, 39, 474-481. Tanaka, K., Horiuchi, S., Sato, T., Zana, N., 1980. The earthquake generating stress in the Western Rift Valley of Africa. J. Phys. Earth., 29: 45-57. Tanaka, K., 1983. Seismicity and Focal mechanism of the volcanic earthquakes in Virunga volcanic region. In: Hamaguchi, H (Editor), Volcanoes Nyiragongo and Nyamuragira, Geophysical Aspects, Tohoku University, Sendai, Japan. Thurber, C.H., 1984. Seismic detection of the summit magma complex of Kilauea volcano, Hawaii. Science 223, 165-167. Tinti, S., Mulargia, F., 1985. Effects of magnitude uncertainties in the Gutenberg Richter frequency magnitude law, Bull. Seism. Soc. Am. 75, 1681-1697. Twesigomwe, E., 1997. Probabilistic seismic hazard assessment of Uganda. Ph.D. thesis, Dept. of Physics, Makerere University, Uganda. Vauchez, A., Barruol, G., Tommasi, A., 1997. Why do continents break-up parallel to ancient Orogenic belts? Terra Nova, 9(2), Oxford, 62-66. Villeneuve, M., 1983. Les sillons tectoniques du Precambrien Superieur dans l?Est du Zaire, comparaison avec les directions du Rift Est-Africain, Bull. Centre Rech. Expl. Prod. Elf-Aquitaine, 7, 163-174. Wafula, M., Zana, N., 1990. Focal mechanism Study of Earthquakes, Revue des 135 Sciences Naturelles, Centre de Rechercherches en Sciences Naturelles, Zaire, Vol.1, No1. Walpersdorf, A., Vigny, C., Ruegg, J.C., Huchon, P., Asfaw, L.M., Kirbash, S.A., 1999. 5 years of GPS observations of the Afar Triple Junction area, Journal Geodynamics 28, 225-36. Wang, G. Z. Q. J., Tim Law, K., 1994. Siting in Earthquake zones, A. A. Balkena, 37-52. Wohlenberg, J., 1967. Seismizit?t der ostafrikanishen Grabenzonen zwischen 4?N and 12?S sowie 23?E und 40?E. Thesis, University of Ludwig-Maximilian, Munchen, 95 pp. Wohlenberg, J., 1968. Remarks on the seismicity of East Africa between 4?N- 12?S and 23?E-40?E. Tectonophysics, 8, 567-577. Yurui, H., Tianzhong, Z., 1997. The K method for estimating earthquake activity parameters and effect of the boundary uncertainty of the source region: Discussion on the zoning method, Earthq. Res. In China, 11, 299-305. Zana, N., 1977. The Seismicity of the Western Rift Valley of Africa and associated problems, Doctorate Thesis, Tohoku University, Sendai, Japan, 177 pp. Zana, N., Hamaguchi, H., 1978. Some characteristics of aftershock sequences in the Western Rift Valley of Africa, Tohoku Geophys. Journ. (Sci. Rep. Tohoku Univ., Ver.5), Geophysics, Vol.25, No 2, p55-72. Zana, N., Tanaka, K., 1981. Focal Mechanism of Major Earthquakes in the Western Rift Valley of Africa. Tohoku Geophys. J., 28, 119-121. Zana, N., Kamba, M., Katsongo, S., Janssen, Th., 1989. Recent seismic activity of 136 the Kivu Province, Western Rift Valley of Africa, Physics of the Earth and Planetary Interiors, 58, 52-60. Zana, N., Horiuchi, S., Murakami, E., Takagi, A., 1990. High frequency earthquakes occurring outside the Western Rift Valley of Africa, Tohoku Geophys. Journ. (Sci. Rep.Tohoku Univ.,Ser.5),Vol.33, No1, pp69- 82. Zana, N., Kavotha, K., Wafula, M., 1992a. Estimate of earthquake risk in Zaire, Tectonophysics, 209, 321-323. Zana, N., Tanaka, K., Kasahara, M., 1992b. Main geophysical features related to the Virunga zone, Western rift, and their volcanological implications Tectonophysics, 209, 255-257. Zana, N., Wafula, M., Lukaya, N., Batabolo, M., 2004. The Kabalo earthquake in D. R. Congo on September 11, 1992: Field observations and damages Quelques r?sultats de Recherches en G?ophysique. Centre de Recherches et Pedagogie appliqu?s ( C. R. P. A), I. P.N, Kinshasa, 77-89. Zhang, H., Thurber, C.H., 2003. Double difference tomography: The method and its application to the Hayward fault, California. Bull. Seis. Soc. Am. 93: 1875-93. 137 APPENDIX A Table A.1: List of earthquakes used in this study for the regression analysis of Mb(USGS) versus M (LWI) Date Origin time Latitude Longitude Depth(km) M(LWI) Mb(USGS) 1965/04/25 10:01:07 2.55S 28.85E 5 4.5 5.05 1965/06/21 11:12:02 4.10S 35.10E 21 6.3 5.05 1966/03/20 01:42:50 0.70N 29.80E 24 7 6.2 1966/03/20 02:39:40 1.10N 30.00E 16 5.4 5.4 1966/03/20 08:55:35 0.80N 29.80E 5 5.6 5.3 1966/03/21 01:30:38 0.80N 29.80E 5 5.3 5.4 1966/03/21 09:23:49 0.80N 29.90E 6 5 5.1 1966/04/01 07:59:41 0.60N 29.50E 17 4.5 4.5 1966/04/06 01:17:54 0.80N 29.90E 41 4.6 4.4 1966/04/07 00:09:11 0.60N 29.90E 33 4.4 4.4 1966/04/14 13:16:19 0.90N 30.00E 33 5 4.8 1966/04/15 03:08:15 0.80N 30.10E 26 5.6 5.2 1966/04/16 14:43:18 0.80N 29.90E 11 5.1 5.3 1966/04/30 20:38:47 0.70N 29.80E 33 4.8 4.9 1966/05/17 07:03:33 0.70N 29.90E 35 5.5 5.3 1966/05/18 01:46:34 0.70N 29.90E 34 4 4.5 1966/05/29 02:26:11 0.70N 30.00E 33 5.2 4.8 1966/06/03 07:14:42 0.60N 29.90E 33 5 5 1966/06/17 18:31:55 0.80N 29.90E 33 5.3 5 1967/08/24 23:14:45 10.50S 27.30E 21 4.8 5 1967/10/30 19:55:38 2.00N 31.30E 33 4.2 5.3 1967/10/30 22:22:39 2.10N 31.20E 33 4.2 4.8 1967/10/31 13:51:04 2.00N 31.20E 33 5.2 5.3 1967/11/11 02:28:46 2.00N 31.50E 33 5.3 5.1 1967/11/13 20:29:05 1.90N 31.70E 33 5.1 5.3 1969/04/22 21:59:11 1.90N 31.50E 28 5 6 1969/04/29 19:54:45 0.80S 30.70E 33 5 5 1969/06/08 15:40:06 6.10S 30.60E 33 4.5 3.8 1969/06/08 16:52:40 6.10S 30.80E 33 5.2 4.7 1969/09/09 16:46:44 2.60S 24.70E 33 5.3 5.2 1970/01/15 08:50:39 11.20S 34.50E 33 5.7 4.9 138 Table A1 (Continued) Date Origin time Latitude Longitude Depth(km) M(LWI) Mb(USGS) 1970/01/19 20:10:20 7.70S 25.90E 37 4.6 4.6 1970/01/27 20:44:45 3.30S 25.80E 55 5.3 4.2 1970/05/10 23:45:22 10.90S 32.70E 33 5.2 3.9 1970/05/23 23:13:59 8.00S 30.80E 33 4.6 4 1970/08/09 00:26:55 6.00S 34.80E 33 5.3 4.8 1971/01/04 15:14:35 3.60N 32.40E 33 4.1 4.4 1971/01/06 04:26:26 1.00S 25.10E 33 4.3 4.4 1971/01/14 20:47:52 11.20S 34.30E 33 4.1 4.4 1971/01/16 09:00:20 1.40S 28.60E 18 5 5 1971/01/27 22:26:16 8.50S 32.20E 33 4.5 3.6 1971/04/18 00:34:34 0.20N 30.10E 4.9 4.6 1971/04/21 18:40:54 0.20N 29.90E 4.5 4.3 1971/11/13 15:47:41 11.00N 39.70E 24 5.6 5.3 1971/11/16 04:27:34 2.00S 26.90E 4.8 4.5 1971/12/31 17:43:11 13.20S 26.40E 4.7 4.6 1972/01/08 17:27:51 0.58N 30.08E 5.2 4.8 1972/02/13 10:02:42 4.50S 34.14E 6.8 5 1972/04/18 15:08:13 0.67N 29.82E 5.2 5 1972/04/25 00:59:51 9.20S 33.30E 6.25 4.7 1972/09/08 04:22:19 9.09S 29.90E 5 4.7 1972/12/12 03:18:45 16.82S 28.06E 6.2 4.7 1972/12/12 03:18:45 16.82S 28.06E 6.2 4.7 1972/12/18 01:18:53 16.71S 28.07E 6.6 5.2 1972/12/27 15:29:38 16.71S 28.05E 5.4 4.6 1973/04/15 13:13:33 7.18S 30.22E 5.9 4.8 1973/04/22 22:03:41 4.10N 31.30E 5.9 4.6 1973/05/03 02:27:15 8.40S 32.60E 5.7 4.5 1973/05/14 20:16:27 6.19S 29.86E 4.6 4.2 1973/05/27 09:26:03 10.53S 34.30E 5.5 4.4 1973/07/07 16:04:11 3.05S 35.60E 5.5 3.8 1973/07/16 18:08:22 10.51S 34.04E 5.4 4.4 1973/08/28 23:43:55 8.76S 31.30E 4.8 4 1973/09/01 11:23:55 13.24S 24.14E 6 4.7 1973/09/21 00:50:30 3.60N 28.00E 4.8 5.1 1973/12/01 16:51:15 0.72N 29.62E 5.4 4.5 1973/12/14 15:04:10 4.71S 28.30E 4.5 4.5 1974/02/10 16:29:26 2.90S 23.30E 4.8 4.7 1974/02/18 09:59:46 3.20S 29.50E 3.9 4.6 1974/04/08 06:12:49 5.50S 35.80E 5.5 4.3 1974/04/24 18:59:53 13.70S 26.60E 5.2 4.4 139 Table A1 (Continued) Date Origin time Latitude Longitude Depth(km) M(LWI) Mb(USGS) 1974/04/25 00:03:49 1.00N 30.10E 5.2 5 1974/08/01 09:36:27 16.70S 28.00E 6.6 5.1 1974/09/23 19:28:17 0.30S 12.90E 7.2 5.9 1974/10/04 17:04:53 11.60S 34.20E 5.3 4.9 1974/10/17 20:04:27 9.50S 32.20E 5.9 4.3 1975/03/26 03:40:48 5.40S 30.20E 5.6 5.1 1975/04/06 04:52:08 5.10S 27.70E 5 4.7 1975/08/06 07:37:30 4.37S 35.90E 6.5 5.1 1976/01/09 08:05:48 5.46S 28.74E 5.3 5.2 1976/02/05 07:46:28 2.97S 36.97E 5.5 4.2 1976/05/15 08:09:57 4.46N 19.35E 7 5.6 1976/07/20 22:56:50 9.00S 28.62E 4.8 4.1 1976/08/23 00:30:11 7.39S 30.87E 4.5 5.2 1977/01/04 20:44:39 7.44S 38.52E 5.7 5.2 1977/01/14 02:52:35 1.71S 28.95E 4 4.8 140 Table A2: List of hypocenters and magnitudes of earthquakes used to determine P receiver functions Date Origin time Latitude Longitude Depth(km) Magnitude Region 2004/07/25 14:35:19 2.427S 103.981E 582 Mb6.8 Sumatra, Indonesia 2004/09/06 12:42:59 55.372S 28.976W 6.9 Mw6.9 South Sandwich island 2004/11/21 11:41:08 15.679N 61.706W 14 Mb6.3 Leeward Islands 2005/01/01 06:25:45 5.099N 92.304E 13 Mw6.7 Sumatra, Indonesia 2005/01/12 08:40:04 0.878S 21.194W 10 Mb5.7 Mid-Atlantic Ridge 2005/01/23 20:10:12 1.198S 119.933E 11 Mb5.8 Sulawesi, Indonesia 2005/02/14 23:38:09 41.728N 79.440E 22 Mw6.1 Xinjiang, China 2005/02/16 20:27:52 36.320S 16.558W 10 Mb6.0 Mid-Atlantic Ridge 2005/02/22 02:25:23 30.754N 56.816E 14 Mb6.0 Central Iran 2005/02/26 12:56:53 2.908N 95.592E 36 Mb6.0 Simeulue, Indonesia 2005/03/14 01:55:56 39.354N 40.890E 5 Mb5.5 Eastern Turkey 2005/03/28 16:09:37 2.085N 97.108E 30 Mw8.6 Sumatra, Indonesia 2005/04/10 10:29:11 1.644S 99.607E 19 Mb6.4 Kepulauan Mentawai 2005/05/14 05:05:18 0.587N 98.459E 34 Mb6.4 Nias, Indonesia 2005/05/19 01:54:52 1.989N 97.041E 30 Mw6.9 Nias, Indonesia 2005/07/05 01:52:03 1.819N 97.082E 21 Mw6.7 Nias, Indonesia 2005/07/24 15:42:06 7.920N 92.190E 16 Mw7.2 Nicobar Islands, India 2005/10/08 03:50:41 34.539N 73.588E 26 Mb6.9 Pakistan 2005/10/29 04:05:56 45.214S 96.898E 8 Mb6.1 SE India Ridge 2005/11/19 14:10:13 2.164N 96.786E 21 Mw6.5 Simeulue, Indonesia 2005/11/27 10:22:19 26.774N 55.858E 10 Mb6.1 Southern Iran 2006/01/02 06:10:50 60.957S 21.606W 13 Mb6.4 South Sandwich island 2006/01/08 11:34:56 36.311N 23.212E 66 Mb6.5 Southern Greece 2006/02/28 07:31:03 28.120N 56.865E 18 Mb5.8 Southern Iran 2006/03/25 07:28:58 27.574N 55.685E 18 Mb5.7 Southern Iran 2006/03/31 01:17:01 33.500N 48.780E 7 Mb5.7 Western Iran 2006/05/16 15:28:26 0.093N 97.050E 12 Mb6.6 Nias, Indonesia 2006/05/26 22:53:59 7.961S 110.446E 13 Mb6.0 Java, Indonesia 2006/06/28 21:02:10 26.925N 55.866E 11 Mb5.8 Southern Iran 2006/07/17 08:19:27 9.284S 107.419E 20 Mb6.1 Java, Indonesia 2006/07/29 00:11:51 37.255N 68.828E 34 Mw5.6 Tajikistan 2006/09/29 13:08:26 10.876N 61.756W 53 Mb5.9 Trinidad and Tobago 2006/11/13 01:26:36 26.052S 63.283W 572 Mb6.3 Santiago, Argentina 141 Date Origin time Latitude Longitude Depth(km) Magnitude Region 2006/12/26 12:26:21 21.799N 120.547E 10 Mb6.4 Taiwan 2007/03/06 05:49:25 0.488S 100.530E 11 Mb5.9 Sumatra, Indonesia 2007/06/02 21:34:58 23.028N 101.052E 5 Mb5.7 Yunnan, China 2007/06/24 00:25:18 55.645S 2.626W 10 Mb5.7 Mid-Atlantic Ridge 2007/08/08 17:05:05 5.859S 107.419E 280 Mb6.5 Java, Indonesia 2007/08/20 22:42:28 8.037N 39.251W 6 Mb6.3 Mid-Atlantic Ridge 2007/09/12 11:10:27 4.438S 101.367E 34 Mb6.9 Sumatra, Indonesia 2007/09/20 08:31:14 1.999S 100.141E 30 Mb6.3 Sumatra, Indonesia 2007/10/24 21:02:51 3.899S 101.020E 21 Mb6.1 Sumatra, Indonesia 2007/11/25 16:02:17 8.277S 118.339E 35 Mb6.2 Sumatra, Indonesia 2007/11/25 19:53:08 8.225S 118.453E 35 Mb6.2 Sumbawa, Indonesia 2007/11/29 19:00:20 14.973N 61.263W 148 Mb6.9 Martinique 142